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IMPLEMENTATION  OF  AN  OPTIMAL  MULTICOMMODITY  NETWORK  FLOW 
ALGORITHM  BASED  ON  GRADIENT  PROJECTION  AND  A  PATH  FLOW  FORMULATION^ 


by 


Dimitri  P.  Bertsckas,  Bob  Goiulron  and  Wei  K.  Tsai* 


ABSTRACT'1' 

The  implementation  of  a  multicommodity  flow  algorithm  into  a  FORTRAN 
code  is  discussed.  The  algorithm  is  based  on  a  gradient  projection  method 
[1]  with  diagonal  scaling  based  on  Hessian  or  Jacobian  information.  The 
flows  carried  by  the  active  paths  of  each  origin-destination  (OD)  pair 
arc  iterated  upon  one  01)  pair  at  a  time.  Active  paths  are  generated  using 
a  shortest  path  algori thm- -one  path  per  OD  pair,  per  iteration.  The  data 
structures  and  memory  retpii rements  of  the  algorithm  are  discussed  and  are 
compared  with  those  of  other  formulations  based  on  link  flows  associated 
with  each  origin,  and  aggregate  link  flows. 
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1.  Optimal  Multicommodity  Flow  Problem  Formulation 


We  have  a  directed  network  with  set  of  nodes  L  and  set  of  links  N. 

Let  W  be  a  collection  of  ordered  node  pairs  referred  to  as  origin-destination 
(OD)  pairs.  For  each  OD  pair  wrW  we  arc  given  a  positive  number  rw  represent¬ 
ing  input  flow  into  the  network  from  origin  to  destination.  Let  be  a 
given  set  of  directed  paths  joining  the  origin  node  and  destination  node 
of  OD  pair  w.  (P  could  be  the  set  of  all  simple  directed  paths  joining 
origin  and  destination,  or  it  could  be  a  restricted  set  of  paths  determined 
a  priori  on  the  basis  of  some  unspecified  considerations).  Note  that  we 
do  not  exclude  the  possibility  that  two  distinct  OD  pairs  have  the  same 
origin  and  destination  and  possibly  a  different  set  of  paths,  but  are  as¬ 
sociated  with  different  classes  or  types  of  traffic. 

Let  Xp  be  the  flow  carried  by  a  generic  path  p.  The  optimization 
variables  of  the  problem  are  x^,  pcP^,  weW  and  must  satisfy  the  constraints 


V  wr.W, 


x  >  0 

P  ~ 


V  pc.P  ,  wr.W . 
1  w 


Let  x  be  the  vector  of  all  path  flows 


x  =  tx  I  peP  ,  wtW 1 
p  1  r  w 


For  each  link  (i,j)  and  OD  pair  w  we  are  given  a  continuously  dif¬ 


ferentiable  function  T^j(x,w),  which  is  to  be  interpreted  as  the  length 
of  link  (i,j)  when  the  path  flow  vector  is  x.  In  data  communication  rout¬ 
ing  and  traffic  assignment  problems  T.j(x,w)  usually  has  the  interpretation  of 
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marginal  delay  and  travel  time  respectively  (see  [1]-[19]).  We  assume 
that  for  all  feasible  x  and  all  weW 

T^Cx.w)  >0  ,  V  (i ,  j)eL  ,  (4) 


The  length  of  a  path  peP^  when  the  path  flow  vector  is  x  is  defined  by 


L  (x.w)  =  l  T  (x,w) 

P  (i.j)ep  J 


(5) 


i.e.  it  is  the  sum  of  lengths  of  its  links. 

Hie  problem  we  are  considering  is  the  following: 

Find  a  path  flow  vector  x*  satisfying  the  constraints  (1),  (2)  and  such 
that  for  every  weW  and  peP^ 


xj  >  0  Lp(x*,w)  <  Lpf(x*,w),  V  p'£Pw. 


(6) 


In  other  words  we  are  looking  for  a  path  flow  pattern  x*  whereby  the 
only  paths  that  carry  positive  flow  are  shortest  paths  with  respect  to  the 
link  lengths  T.j(x*,w). 

Hie  problem  described  above  includes,  among  others,  problems  of 
optimal  routing  in  data  networks  [l]-[8]  and  (possibly  asymmetric)  traffic 
assignment  problems  in  transportation  networks  [9] -[19].  We  refer  to  the 
references  just  cited  for  extensive  discussions.  The  survey  paper  [1] 
describes  in  detail  the  data  communication  context.  A  typical  formulation 
there  is  to  find  a  feasible  path  flow  vector  x  that  minimizes 


l 

(i  »j) 


D. .(F.  .) 
ij  ir 


(7) 


where  is  a  monotonically  increasing,  twice  differentiable  function  of 
the  total  flow  of  the  link  (i,j)  given  by 


Fii  =  l  l  x  <S(P.i.j) 
J  weW  peP  1 


(8) 


where 


6  (p ,  i ,  j  }  = 


1  if  link  (i,j)  belong  to  path  p 


0  otherwise. 


(9) 


It  can  be  shown  (see  e.g.  [1])  that  if  we  make  the  identification 


T. .  =  D! .  :  The  first  derivative  of  D. . 

1J 


(10) 


the  routing  optimization  problem  falls  within  the  framework  of  the  general 
multicommodity  flow  problem  described  earlier. 
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A  Projection  Method  for  Solving  the  Multicommodity  Flow  Problem 
The  MULTIFLO  and  MULTIFL01  codes  given  in  Appendices  I  and  II  of  this 
report  implement  an  algorithm  that  solves  the  problem  of  the  previous 
section  for  the  case  where  for  all  OD  pairs  weW 

Pw  =  Set  of  all  simple  paths  joining  the  origin  and  destination  of  w. 

The  set  of  OD  pairs  is  divided  into  C  groups  called  commodities.  All  OD 
pairs  of  a  commodity  have  the  same  origin  node.  Furthermore  the  data  structures 
of  the  codes  can  handle  only  the  case  where  the  lengths  T „  (x,w)  depend  on 
w  through  the  corresponding  commodity  c.  That  is 


Tjj(x,w)  =  T^Cx.w),  V  (i,j)fl-,  and  OD  pairs  w,  w  of  the  same 

commodity  c. 

It  is  also  assumed  that  for  all  feasible  F 
3T.  . 

— >0  V  (i,j)  belonging  to  the  path  p 
XP 

MULTIFLO  and  MULTIFLOl  operate  as  follows: 

At  the  beginning  of  the  kth  iteration  wc  have  for  the  generic  OD 
pair  w  W  a  set  of  active  paths  Pw  consisting  of  at  most  (k-1)  paths. 
(These  paths  were  generated  in  earlier  iterations  and  it  is  implicitly 
assumed  that  all  other  paths  carry  zero  flow).  The  following  calculation 
is  executed  sequentially  for  each  commodity--first  for  commodity  1,  then 
for  commodity  2,  and  so  on  up  to  the  last  commodity  C: 
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Step  1 :  A  shortest  path  that  joins  the  origin  node  for  the  commodity 
with  all  other  nodes  is  calculated.  The  length  for  each  link  (i,j)  used 
for  this  calculation  is  T  (x,w)  where  x  is  the  current  path  flow  vector. 
These  shortest  paths  are  added  to  the  corresponding  list  of  active  paths 
of  each  OD  pair  of  the  commodity  if  they  are  not  already  there,  so  now 
the  list  of  active  paths  for  each  OD  pair  of  the  commodity  contains  at 
most  k  paths. 

Step  2:  Each  OD  pair  w  of  the  commodity  is  taken  up  sequentially.  For  each 
active  path  p  of  w  the  length  L  Icf.  (5)]  is  calculated  together  with  an 
additional  number  a  called  the  stepsizc  (more  on  the  choice  of  this 

p  — i - 

later) .  Both  L  and  a  are  calculated  on  the  basis  of  the  current  total 
P  P 

link  flow  vector.  Let  p  be  the  shortest  path  calculated  in  Step  1  for 
the  OD  pair.  The  path  flows  of  all  paths  p  ^  p  are  updated  according  to 


max  {0,  x  -a  (L  -L-) }  if  L  >  L~ 
P  P  P  p  P  P 


otherwise. 


The  path  flow  of  the  shortest  path  p  is  then  adjusted  so  that  the  sum 
of  flows  of  all  active  paths  equals  r^  as  required  by  the  constraint  (1), 


x—  «-  r  -  )  x  . 

p  w  ‘  .—  p 

1  active  pfp  1 


In  other  words  an  amount  x  or  a  (L  -L— )  is  shifted  from  each  nonshortest 

P  p  p  p 

path  to  the  shortest  path  p--whichever  is  smaller.  The  total  link  flows  F. 

are  adjusted  to  reflect  the  changes  in  x  and  x— . 

P  P 

The  rationale  for  iteration  (11)  is  explained  in  [1],  {6J,  [8],  [9]. 
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It  is  based  on  a  gradient  projection  method  [9],  [21].  Note  that  it  is 

possible  that  <  L—  for  some  p  /  p  even  though  p  was  calculated  earlier 

as  a  shortest  path.  The  reason  is  that  by  the  time  L  and  L—  are  computed 

P  P 

the  total  link  flow  vector  may  have  changed  since  the  time  the  shortest 
path  has  been  calculated  due  to  iterations  on  the  path  flows  of  other  OD 
pairs  of  the  same  commodity. 

Regarding  the  choice  of  the  stepsize  a  ,  the  MULTIFLO  and  MULTIFLOl 
codes  use  the  following  formula  for  all  p  ^  p 


a 


P 


where 


S 

P 


(i 


l 

,j)cLp3 


and  L  is  the  set  of  links 
P 


(13) 


(14) 


=  ((i.j)j  (i, j)  belongs  to  either  p  or  p  , 
but  not  to  both  p  and  p7 . 

The  rationale  for  this  is  as  follows: 

If  we  interpret  the  algorithm  as  one  that  tries  to  satisfy  the  equation 

(16) 

L  -  L—  =  0,  Vpwithx  >0, 

P  P  P 


a  natural  choice  for  a  is 

P 


a 

P 


A(L 


V 


(17) 


where  A(L  -L— )  is  the  variation  of  (1,  -L— )  resulting  from  a  small  variation 
P  P  P  P 


(i 


«')x 

P 


thereby  justifying  the  use  of  the  stepsize  (13),  (14). 

If  one  wishes  to  employ  the  formula  (IS)  for  the  stepsize  it  is 
necessary  to  modify  the  codes.  These  modifications  should  not  be  too 
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difficult  for  an  experienced  user.  Another  possibility  is  to  use  a  smalier 
value  of  stepsize  ctp  than  the  one  given  by  (13)--for  example  ot^  =  pS^ 
pe(0,l)  is  a  fixed  relaxation  parameter.  (A  smaller  stepsize  enhances  the 
convergence  properties  of  the  algorithm  but  may  deteriorate  its  rate  of 
convergence).  This  can  be  accomplished  without  any  changes  in  the  code 

by  simply  introducing  the  relaxation  parameter  p  in  the  subroutine  that 

3T.  . 

calculates  — [cf.  (14)]. 

3xp 

In  the  MULTIFLO  code  a  shortest  path  tree  is  generated  and  stored 

at  each  iteration  for  each  commodity.  As  a  result  the  memory  storage  for 

shortest  paths  is  proportional  to  the  number  of  iterations  so  for  large 

problems  one  cannot  execute  a  large  number  of  iterations  without  incurring 

a  heavy  penalty  for  disk  I/O.  MULTIFLO  will  usually  find  in  five  to  ton 

iterations  what  is  for  most  practical  problems  an  adequate  approximation 

to  an  optimal  solution.  This  is  particularly  true  of  lightly  loaded  net¬ 
works  (e.g.  with  utilization  of  all  links  less  than  60",,  at  the  optimum). 

For  heavily  loaded  networks  the  number  of  required  iterations  usually 

tends  to  be  larger  (say  10-30).  It  should  be  a  rare  occasion  when  a  user 

will  require  more  than  thirty  iterations  for  his  practical  problem. 

MULTI FL01  differs  from  MULTIFLO  only  in  the  method  used  for  storing 
the  active  paths.  MULTIFL01  stores  explicitly  all  active  paths  in  a 
single  array  rather  than  storing  them  implicitly  through  the  generated 
shortest  path  trees.  As  a  result  the  memory  storage  of  MULTIFL01  depends 
on  the  number  of  active  paths  generated  and  is  largely  independent  of 
the  number  of  iterations  executed.  For  certain  problems  including 
situations  where  a  large  number  of  iterations  is  desired  MULTIFL01  may 
hold  astorage  advantage  over  MULTIFLO.  Both  codes  generate  identical 
numerical  results  although  MULTIFL01  appears  to  be  somewhat  faster  on 
sample  test  problems. 
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3.  Data  Structures  for  Representing  the  Problem 

The  data  structures  of  MULTIFLO  and  MULTIFL01  are  described  in  the 
code  documentation.  The  problem  input  structure  will  be  illustrated  here 
by  means  of  the  5  node-6  link  network  shown  in  Figure  l: 


•o 


These  arrays  specify  the  network  topology. 


FRSTOU (NODE) :  The  first  link  out  of  NODE 
LASTOU(NODE) :  The  last  link  out  of  NODE 
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NODE 

FRSTOU 

LASTOU 

1 

1 

1 

2 

2 

4 

5 

3 

3 

3 

4 

6 

6 

S 

0 

0 

Note  that  all  arcs  with  the  same  head  node  must  be  grouped  together  in 
the  arc  list.  A  node  with  no  outgoing  links  is  recognized  via  FRSTOU  =  0 

Arc  Length  Arrays  (STARTNODE,  ENDNODE) 

These  arrays  also  specify  the  network  topology: 

STARTNODE  (ARC) :  The  head  node  of  ARC 
ENDNODE  (ARC):  The  tail  node  of  ARC 


ARC 

STARTNODE 

ENDNODE 

1 

1 

3 

1 

4 

3 

3 

5 

4 

2 

3 

5 

2 

4 

6 

4 

5 

Commodity  Length  Arrays  (ORGID,  STARTOD) 

ORGID  (COMMODITY):  The  origin  node  of  COMMODITY 

STARTOD  (COKWODITY) :  A  pointer  to  the  first  OD  pair  of  COMMODITY  on 


the  OD  pair  list 


T 
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For  the  example  of  Figure  1  we  will  assume  three  commodities 


COMMODITY 

ORGID 

STARTOD 

1 

2 

1 

2 

1 

3 

3 

1 

4 

Note  that  it  is  required  that  OP  pairs  are  listed  sequential ly  by  com¬ 
modity,  i.e.  the  OD  pairs  of  commodity  1  are  listed  first,  followed  by 
the  OD  pairs  of  commodity  2,  etc.  Therefore  the  STARTOt)  array  together 
with  the  total  number  of  OD  pairs  specify  all  OD  pairs  associated  with 
each  commodity. 

OD  Pair  Length  Arrays  (PEST,  INPUT, FLOW) 

DEST(OD):  The  destination  node  of  OD 
INPUTFLOW(OD) :  The  input  traffic  of  OD 


OD  DEST  INPUTJ^LOW 


1 

3 

problem  dependent 

2 

5 

It 

3 

3 

It 

4 

4 

tl 

S 

5 

ft 

From  the  arrays  ORGID,  STARTOD  and  DEST  together  with  the  total  number 
of  OD  pairs  the  set  of  OD  pairs  corresponding  to  each  commodity  is  com¬ 
pletely  specified.  For  our  example  these  arc: 


i 
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COMMODITY 

OD  PAIRS 

1 

(2,3),  (2,5) 

2 

(1,3) 

3 

(1,4),  (1,5) 

Additional  input  information  is  required  to  calculate  the  link 

9T.  . 

lengths  T.  .  and  their  first  derivatives  in  the  subroutine  DERIVS 

1J  XP 

and  DERIV1 .  This  is  of  course  problem  dependent.  The  listing  of  Appendix  I 
gives  an  example  which  is  typical  of  routing  problems  in  data  networks  [cf. 
equations  (7) -(10)]. 
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4 .  Memory  Requirements  -  Comparisons  with  Other  Methods 

The  memory  storage  requirements  of  both  MULTIFLO  and  MULTIFL01  are 
substantial,  but  this  is  true  for  all  methods  that  provide  as  output 
not  only  the  optimal  total  link  flows  but  also  detailed  information 
about  the  optimal  routing  from  origins  to  destinations  (i.e.  optimal 
path  flows) . 

Assuming  that  1  byte  is  allocated  for  a  logical  variable,  2  bytes 
are  allocated  for  storing  a  node  or  link  identification  number  and  an 
iteration  number,  4  bytes  are  allocated  for  storing  a  commodity,  OD  pair 
or  path  identification  number,  and  4  bytes  are  allocated  for  storing  a 
real  number  (e.g.  a  path  or  link  flow)  the  total  array  storage  in  bytes 
of  MULTIFLO  during  execution  is 

6nN  +  9n^  +  6nc  +  6ngD  +  10np  +  2nj  nN  nc  (19) 


where : 

n„:  Number  of  nodes 
N 

n^:  Number  of  links 

n^:  Number  of  commodities 

n^:  Number  of  OD  pairs 

np:  NurabeT  of  active  paths  generated 

n^:  Number  of  iterations. 

Additional  storage  is  required  for  information  necessary  to  calculate 
link  lengths  and  their  derivatives  but  this  is  typically  of  order  O(n^) 
and  is  not  significant. 

The  dominant  array  as  far  as  storage  of  MULTIFLO  is  concerned  is  the 
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triple  indexed  PRED  array  which  stores  the  shortest  path  trees  generated 
for  each  commodity  at  each  iteration.  This  array  accounts  for  the  last 
term  2njn^n^,  in  (19).  The  term  10np  is  also  substantial  since  the 
number  of  active  paths  np  can  be  as  large  as  n^gg.  However,  because 
the  algorithm  stores  a  path  only  once  at  the  iteration  it  is  first  gen¬ 
erated  and  does  not  duplicate  it  if  it  is  generated  again  later,  the 
actual  number  np  is  typically  much  smaller  than  n^ngg.  This  was  con¬ 
firmed  by  extensive  computational  experimentation,  that  showed  that  except 
for  very  heavily  loaded  networks  the  actual  number  of  active  paths  Op  was 
typically  no  more  than  2ngg(!)  and  often  considerably  less.  We  conclude 
therefore  that  the  dominant  bottleneck  for  storage  is  the  shortest  path 
description  array  PRED  requiring  2n jn^n^  bytes. 

In  the  MULTIFL01  code  the  array  PRED  is  not  used.  In  its  place  the 

array  PDESCR  is  used  which  requires  storage  of  2npnN  at  most.  This 

calculation  assumes  conservatively  that  a  path  has  n^  links.  However 

in  practice  the  actual  storage  for  PDESCR  is  several  times  less  than 

2npn^.  If  we  adopt  the  rough  estimate  np  -  2ngp  then  we  conclude  that 

the  storage  requirements  of  MULTIFLO  and  MULTIFL01  are  roughly  comparable 

n0D 

if  the  number  of  iterations  n^  is  comparable  to  something  between  - — 

n0D  n0D  C 

and  ^ —  with  MULTIFL01  becoming  definitely  preferable  if  nj  =  - — . 

C  C 

MULTIFL01  is  also  preferable  for  problems  that  are  solved  repetitively  with 

minor  variations  in  their  data  since  then  the  knowledge  of  the  path 
description  array  PDESCR  can  be  fruitfully  exploited.  This  is  not  pos¬ 
sible  with  MULTIFLO. 

In  large  problems  where  only  the  total  link  flows  are  of  interest 
(e.g.  traffic  assignment  problems)  a  different  algorithm  [e.g.  the  flow 
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Deviation  ior  the  Frank-Wolfe)  method  [3],  [8]  or  the  Cantor-Gerla  (or 
simplicial  approximation)  method  [4],  [18],  may  he  preferable  over 
MULTIFLO  or  MULTIFL01,  since  then  storage  of  order  0(n  )  or  perhaps 

L 

0(n j,n^)  is  required.  However  when  detailed  routing  information  is  of 
interest  the  memory  storage  requirements  of  MULTIFLO  are  competitive 
with  those  of  other  methods  based  on  shortest  paths  including  the  Flow 
Deviation  and  Cantor-Gerla  methods.  The  reason  is  that  detailed  rout¬ 
ing  information  can  be  provided  by  these  methods  only  if  the  shortest 
paths  generated  at  each  iteration  are  stored  explicitly  in  an  array 
such  as  PRED,  and  as  mentioned  earlier  this  is  the  main  memory  storage 
bottleneck . 

There  are  algorithms  that  can  solve  multicommodity  flow  problems 
and  provide  detailed  routing  information  without  requiring  the  generation 
and  storage  of  shortest  paths.  These  algorithms  are  based  on  a  link  flow 
formulation  [20],  or  the  link  flow  f rac t i on  formulation  due  to  Gallager 
[2],  [5],  [7]  whereby  the  optimization  variables  are  the  flows  or  fractions 
of  flow  respectively  for  each  commodity  that  arc  routed  along  each  link. 

Tbe  storage  requirement  for  these  algorithms  is  of  order  0(n  n  )  and  is 
independent  of  the  number  of  iterations.  Whei  we  compare  this  storage 
with  the  Ofn^n^n^)  storage  of  algorithms  based  on  shortest  paths  we  see 
that  link  flow  formulations  hold  an  advantage  in  terms  of  storage  for 
problems  where  a  large  number  of  iterations  is  desirable.  The  reverse 
is  true  if  the  number  of  iterations  required  for  adequate  solution  of 
the  problem  is  small,  or  if  the  number  of  links  is  much  larger  than  the 
number  of  nodes. 


We  finally  note  a  final  advantage  of  the  path  flow  formulation  over  link 
flow  formulations.  When  the  set  of  paths  for  each  OD  pair  is  restricted  to 
be  a  given  strict  subset  of  the  set  of  all  possible  simple  paths  it  is  extremely 
cumbersome  to  use  a  link  flow  formulation.  By  contrast  it  is  straightforward 
to  modify  the  MULTIFL01  code  to  handle  this  situation. 


i 
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APPENDIX  I:  MULTIFLO  Code 

The  following  FORTRAN  code  works  on  the  VAX  family  of  computers.  It 
consists  of  a  DRIVER  program  and  several  subroutines: 

LOAD:  Reads  network  topology  and  link  length  data  from  disk. 

MULTIFLO:  This  is  the  main  algorithm. 

SP:  Calculates  a  shortest  path  tree  from  an  origin  node  to  all  other 
nodes. 


PRFLOW:  Prints  out  to  disk  problem  data  and  algorithmic  results. 

DERIVS:  This  user  supplied  routine  calculates  for  a  given  link  (i,j)  its 

9T.  . 

length  T  (D1CAL)  and  the  length  derivative  — (D2CAL) . 
ij  3  P 

DERIV1 :  This  routine  is  the  same  as  DERIVS  except  that  it  calculates 

3T.  . 

the  length  T.^  (D1CAL)  but  not  the  length  derivative  . 

DELAY:  This  user  supplied  routine  is  useful  only  if  the  multicommodity 

flow  problem  is  a  routing  optimization  problem  of  the  form  (7) - (10) 
as  described  in  Section  1.  For  asymmetric  traffic  assignment  problems 
it  has  no  purpose.  It  calculates  the  total  delay 


y  D .  .  (F .  . ) 

,  -  L  u  ij 


where  D!  .  =  T.  .  [cf.  (7) -(10)1.  The  value  of  D.  . (F..)  is  calculated 
l  j  l  j  1  ’ 1  ij  v  ij 

using  the  function  DCAL. 

Two  versions  of  the  shortest  path  routine  SP  are  provided  (SHORTPAPF. 
and  SHORTHEAP)  which  can  be  used  interchangeably.  SHORTHEAP  is  recommended 
for  problems  where  there  are  only  few  destinations  for  each  commodity. 
Otherwise  SHORTPAPE  based  on  [23]  should  be  preferable. 

A  program  (SETUP)  is  also  provided  for  the  purpose  of  creating  the  data 
describing  the  problem  in  a  format  that  is  compatible  with  the  LOAD  routine- 
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The  routines  LOAD,  DERIV1 ,  DERIVS,  DELAY,  and  DCAL  supplied  in  this 
appendix  correspond  to  the  most  commonly  solved  optimal  routing  problem  in 
data  communication  network  applications  whereby  a  capacity  C.^  is  given  for 
each  link  (i,j)  (this  is  the  array  BITRATE  in  the  code)  and 

(M/M/1  Queueing  Delay)  (A . 1 ) 


D.  .(F.  .) 
ij  iJ 


F.  . 

IL 


C.  ,-F.  . 
ij  iJ 


T.  . (F. . )  = 

ij  ij 


3T. .(F.  .) 
ij  v  ir 

3f.  . 

ij 


C.  . 


(C. .-F.  .)‘ 
ij  i y 


2C.  . 

_ LL 


(C.  .-F.  .)' 
ij  ij 


Because  D^(F^^)-*-oo  as  F. ^  ->  C.  .  these  formulas  have  been  modified  so  that 
if  Cjj  »  where  pc  (0,1)  is  a  parameter  set  by  the  user,  then  D_  , 

3T\. 

T--,  gp1^-  are  calculated  using  a  quadratic  function  which  has  the  same 
J  ij  F.  . 

value,  first  and  second  derivatives  as  - — ~ —  at  the  breakpoint  pCL  . . 

ij'  ij  1J 

In  the  program  the  parameter  p  is  given  by  the  variable  MAXUTI  set  in 

the  subroutine  LOAD  to  0.99.  The  user  may  wish  to  change  this  value. 

The  guideline  is  that  p  should  be  set  at  a  value  exceeding  the  maximum 

link  utilization 

F.  . 

max  -TT-1- 
(i.j)eL  ij 


at  the  optimal  solution.  This  trick  gets  around  situations  whereby  the 
input  flows  are  so  large  that  exceeding  some  of  the  link  capacities  dur¬ 
ing  some  phase  of  the  algorithm  is  inevitable. 


-J1  - 


The  MULTI  FLO  code  will  stop  computing,  when  one  of  two  conditions  is 
met:  Either  the  maximum  number  of  iterations  (MAX ITER)  is  exceeded  or  a 
normalized  measure  of  deviation  from  the  optimal  solution  falls  below  a 
certain  tolerance  (T0L1 .  litis  measure  is  roughly  equal  to  the  percentage 
of  input  traffic  of  an  01)  pair  that  does  not  lie  on  a  shortest  path  (maxi 
mi  zed  over  all  OD  pairs),  and  its  magnitude  is  not  substantially  affected 
by  the  size  of  the  problem,  both  convergence  parameters  MAXITER  and  TOL 
are  set  by  the  user  in  the  subroutine  LOAD. 


/ 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C  DRIVER 

c 

C  'DRIVER'  IS  A  SIMPLE  EXECUTIVE  TO  INVOKE  THE  'MULTIFLO'  COMMODITY 

C  ROUTING  PROGRAM.  'DRIVER*  INVOKES  SUBPROGRAM  'LOAD'  TO  READ 

C  DATA  INTO  'MULTIFLO'  INPUT  COMMON  BLOCKS.  FILES  READ  BY 

C  'LOAD'  ARE  CREATED  BY  A  TERMINAL  SESSION  WITH  THE  USER  FOR 

C  NETWORK  DEFINITION  THROUGH  THE  USE  OF  PROGRAM  'SETUP'. 

C 

C  EXECUTION  STEPS  FOR  PROGRAM  'DRIVER' 

C 

C  1)  ASSIGN  FORTRAN  UNIT  01  AS  CREATED  BY  PROGRAM  'LOAD* 

C  2)  ASSIGN  FORTRAN  UNIT  02  AS  CREATED  BY  PROGRAM  'LOAD' 

C  3)  ASSIGN  FORTRAN  UNIT  06  AS  A  DESIGNATED  OUTPUT  FILE 

C 

C  E .G. : 

C  $  ASSIGN  NETWORK . DAT  FOROOl 

C  $  ASSIGN  TRAFFIC.DAT  FOR002 

C  §  ASSIGN  OUTPUT . DAT  FOR006 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

PROGRAM  DRIVER 

LOAD  FORTRAN  UNIT  01  AND  FORTRAN  UNIT  02  FROM  DISK  AS  CREATED 
FROM  PROGRAM  'SETUP' 

INCLUDE  'PARAM.DIM' 

INCLUDE  ' PATHS . BLK ' 

INCLUDE  'NETWRK.PRM' 

INCLUDE  ' CONVRG . PRM ' 

INTEGER  COMMOD I TY. ORIGIN. DE STOD , OD , PATH 
CALL  LOAD 

EXECUTE  THE  'MULTIFLO'  NETWORK  ALGORITHM.  'MULTIFLO'  SCHEDULES 
ITS  OWN  OUTPUTS  TO  FORTRAN  UNIT  06  ON  EACH  ITERATION 

INITIALIZE  THE  TIMER 
CALL  LIB$INIT_TIMER 
CALL  MULTIFLO 

RECORD  THE  COMPUTATION  TIME 
CALL  LIB$SHOW_TIMER 

PRINT  MAX  LINK  UTILIZATION  (RELEVANT  FOR  M/M/1  QUEUEING  DELAY 
OPTIMIZATION) 

UMAX=0 .0 
DO  100  1=1, NA 

UMAX=MAX (UMAX,FA(I)/BITRATE(I)) 

00  CONTINUE 

WRITE (6,*) 'MAXIMUM  LINK  UTILIZATION’ 

WRITE (6,*)  UMAX 

PRINT  FINAL  PATH  FLOW  INFO 

WRITE  (6,  *)  'ORIGIN  /  DESTINATION  /  PATH  #  /  PATH_FLOW ’ 

DO  1000  COMMODI TY= 1 , NUMCOMMOD 
ORIGIN=ORGID (COMMODITY) 

DO  500  OD=STARTOD (COMMODITY) , STAR TOD (COMMODITY* 1 ) - 1  _ 


2. 


DESTOD=DEST (OD) 

PATH=OD 

DO  WHILE  (PATH.GT.O) 

WRITE (6, *) ORIGIN, DESTOD , PATH , FP (PATH) 
PATH=NEXTPATH (PATH) 

END  DO 
CONTINUE 
CONTINUE 
STOP 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C  LOAD 

c 

C  'LOAD'  READS  IN  DATA  FROM  DISK  CREATED  WITH  PROGRAM  ’SETUP’  FOR 

C  USE  BY  PROGRAM  ’MULTIFLO’.  NETWORK  SPECIFICATION  DATA  RESIDES 

C  ON  FORTRAN  UNIT  01  AND  NETWORK  TRAFFIC  SPECIFICATION  DATA 

C  RESIDES  ON  FORTRAN  UNIT  02. 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

SUBROUTINE  LOAD 
IMPLICIT  NONE 


C 

C 


********************* 


INCLUDE  COMMON  BLOCKS  ********************* 


C 


c 

c. 

c 

c 


INCLUDE  ’ PARAM .DIM ’ 

INCLUDE  ’ NETWRK . PRM * 

INCLUDE  ' CONVRG . PRM ’ 

********************  LOCAL  VARIABLE  DEFINITIONS  ****************** 

INTEGER  I 

DO  LOOP  INDEX 


C 

C  *********************  EXECUTABLE  code  *************************** 

C  TERMINATION  PARAMETERS.  MAXITER  GIVES  THE  MAX  #  OF  ITERATIONS 

C  TOL  IS  A  SOLUTION  ACCURACY  TOLERANCE.  RECOMMENDED  VALUES 

C  ARE  0.01  TO  0.0001.  THE  PROPER  VALUE  OF  TOL  IS  LARGELY 

C  INDEPENDENT  OF  THE  PROBLEM  SIZE. 

MAXITER=20 
TOL=0 . 0 1 

C  THE  FOLLOWING  PARAMETER  MAKES  SENSE  ONLY  FOR  ROUTING  PROBLEMS 

C  WHERE  AN  M/M/1  QUEUING  FORMULA  IS  USED  FOR  DELAY. 

C  IT  GIVES  THE  THRESHOLD  FRACTION  OF  CAPACITY  BEYOND  WHICH 

C  THE  DELAY  FORMULA  IS  TAKEN  TO  BE  QUADRATIC. 

MAXUTI=0 . 99 
C 

C  LOAD  THE  NETWORK  CONFIGURATION  FROM  FORTRAN  UNIT  01 

C 

C  NODE  SPECIFICATIONS 

C 

READ (1 , * ) NN 
DO  1=1. NN 

READ ( 1 . * ) FRSTOU ( I ) .LASTOU(I) 

END  DO 
C 

C  LINK  SPECIFICATIONS 

C 

READ{1,  *)NA 
C 

C  BITRATE  (I )  IS  A  PARAMETER  ASSOCIATED  WITH  LINK  I.  IN  THE 

C  DATA  NETWORK  ROUTING  CONTEXT  IT  HAS  THE  MEANING  OF 

C  TRANSMISSION  CAPACITY  OF  LINK  I. 

C 

DO  1=1, NA 

READ  ( 1 , * ) STARTNODE ( I ) , ENDNODE ( I ) .BITRATE (I) 

END  DO 


C 

C 


INPUT  COMMODITY  DATA  FROM  FORTRAN  UNIT  02 


4.  * 


READ  ( 2 ,  * )  NUV.COMMOD 
DO  I=1,NUMCOXMOD 

READ  (2 ,  * )  ORGID  (I )  ,STARTOD(I) 
END  DO 

READ ( 2 , * ) NUMODPAIR 
DO  1=1 , NUMODPAIR 

READ (2 ,  * ) DEST (I)  , INPUT_FLOW (I) 
END  DO 
RETURN 
END 


n  n  o 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C  MULTI FLO 

c 

C  MULTICOMMODITY  FLOW  ALGORITHM  BASED  ON  A  PATH  FLOW  FORMULATION 

C  UPDATES  THE  PATH  FLOWS  OF  OD  PAIRS  ONE  AT  A  TIME  ACCORDING  TO 

C  AN  ITERATION  OF  THE  PROJECTION  TYPE. 

C 

C  DEVELOPED  BY  DIMITRI  BERTSEKAS,  BOB  GENDRON,  AND  WEI  K  TSAI 

C 

C  BASED  ON  THE  PAPERS: 

C 

C  1)  BERTSEKAS, D. P. ,  "A  CLASS  OF  OPTIMAL  ROUTING  ALGORITHMS 

C  FOR  COMMUNICATION  NETWORKS",  PROC.  OF  5TH  I TERNAT I ONAL 

C  CONFERENCE  ON  COMPUTER  COMMUNICATION  (ICCC-80) , 

C  ATLANTA,  CA.,  OCT.  1980,  PP. 71-76. 

C 

C  2)  BERTSEKAS, D.P.  AND  GAFNI ,E.M. ,  "PROJECTION  METHODS 

C  FOR  VARIATIONAL  INEQUALITIES  WITH  APPLICATION  TO 

C  THE  TRAFFIC  ASSIGNMENT  PROBLEM",  MATH.  PROGR.  STUDY, 17, 

C  D.C. SORENSEN  AND  J.-B.  WETS  (EDS),  NORTH-HOLLAND , 

C  AMSTERDAM, 1982,  PP .  1 39-159. 

C 

C  3)  BERTSEKAS, D.P.,  "OPTIMAL  ROUTING  AND  FLOW  CONTROL 

C  METHODS  FOR  COMMUNICATION  NETWORKS",  IN  ANALYSIS  AND 

C  OPTIMIZATION  OF  SYSTEMS,  (PROC.  OF  5TH  INTERNATIONAL 

C  CONFERENCE  ON  ANALYSIS  AND  OPTIMIZATION,  VERSAILLES, 

C  FRANCE),  A.  BENSOUSSAN  AND  J.L.  LIONS  (EDS), 

C  SPRINGER-VERLAG,  BERLIN  &  NY, 1982,  PP.  615-643. 

C 

C  4)  BERTSEKAS, D.P.  AND  GAFNI ,  E.M.,  "PROJECTED  NEWTON 

C  METHODS  AND  OPTIMIZATION  OF  MULTICOMMODITY  FLOWS", 

C  IEEE  TRANSACTIONS  ON  AUTOMATIC  CONTROL,  DEC.  1983. 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

C 

SUBROUTINE  MULTIFLO 
C 

IMPLICIT  NONE 


*************** 


INCLUDE  COMMON  BLOCKS 


**************************** 


INCLUDE  'PARAM.DIM' 

INCLUDE  ' NETWRK . PRM ' 

INCLUDE  ' CONVRG . PRM  * 

INCLUDE  ' PATHS. BLK' 

C 

C  NODE  ARRAYS  (LENGTH  NN) : 

C 

C  FRSTOU (NODE)  -  FIRST  ARC  OUT  OF  NODE 

C  LASTOU (NODE)  -  LAST  ARC  OUT  OF  NODE 

C  NOTE:  THE  ARC  LIST  MUST  BE  ORDERED  IN  SEQUENCE  SO 

C  THAT  ALL  ARCS  OUT  OF  ANY  NODE  ARE  GROUPED  TOGETHER 

C 

C  ARC  ARRAYS  (LENGTH  NA)  : 

C 

C  FA (ARC)  -  THE  TOTAL  FLOW  OF  ARC 

C  STARTNODE (ARC)  -  THE  HEAD  NODE  OF  ARC 

C  ENDNODE (ARC)  -  THE  TAIL  NODE  OF  ARC 

C 


no  n n n n o n n n n n n n n n n o o o n n o o o o n o o n n 


COMMODITY  LENGTH  ARRAYS  (LENGTH  NUMCOMMOD) : 

ORGID (COMMODITY)  -  THE  NODE  ID  OF  THE  ORIGIN  OF  COMMODITY 
STARTOD (COMMODITY)  -  THE  STARTING  OD  PAIR  IN  THE  ODPAIR  LIST 
CORRESPONDING  TO  THE  ORIGIN  IN  POSITION  RANK 
NOTE:  THIS  SCHEME  ASSUMES  THAT  OD  PAIRS  ARE  LISTED  IN  SEQUENCE 
I.E.  THE  OD  PAIRS  CORRESPONDING  TO  THE  COMMODITY  ONE 
ARE  LISTED  FIRST.  THEY  ARE 

FOLLOWED  BY  THE  OD  PAIRS  OF  THE  COMMODITY  TWO 
AND  SO  ON. 


C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 

C 


ODPAIR  ARRAYS  (LENGTH  NUMOD) : 

DEST (OD)  -  GIVES  THE  DESTINATION  OF  ODPAIR  OD 
INPUT_FLOW (OD)  -  GIVES  THE  INPUT  TRAFFIC  OF  ODPAIR  OD 


PATH  ARRAYS  (LENGTH  DYNAMICALLY  UPDATED) : 

PAIHID (PATH)  -  THE  ITERATION  #  AT  WHICH  PATH  WAS  GENERATED 
NEXTPATH (PATH)  -  THE  NEXT  PATH  FOR  THE  SAME  OD  PAIR  FOLLOWING 
PATH.  IT  EQUALS  0  IF  PATH  IS  THE  LAST  FOR  THAT  OD  PAIR 
FP (PATH)  -  THE  FLOW  CARRIED  BY  PATH 


PATH  DESCRIPTION  LIST  ARRAY  (LENGTH  MAX I TER  *  NUMCOMD * NN) 

PRED (NODE. I TER. COMMODITY)  -  THIS  TRIPLE  INDEXED  ARRAY  SPECIFIES  THE 
SHORTEST  PATH  TREE  GENERATED  AT  ITERATION  ITER 
&  CORRESPONDING  TO  THE  ORIGIN  ASSOCIATED  W/  COMMODITY 
IT  GIVES  THE  LAST  ARC  ON  THE  SHORTEST  PATH  FROM  ORIGIN  TO  NODE. 


*************** 


LOCAL  VARIABLE  DEFINITIONS  ************************ 


INTEGER*  2  PRED (NNN.NMAXITER.NNORIG) 

PATH  DESCRIPTION  ARRAY  -  CONTAINS  SHORTEST 
PATH  TREES  FOR  ALL  ITERATIONS 

LOGICAL  SPNEW 

LOGICAL  INDICATING  A  NEW  PATH  FOUND 
LOGICAL  SAME 

LOGICAL  INDICATING  A  NEW  SHORTEST  PATH  AT .READY  EXISTING 
INTEGER  NODE 

NODE  IDENTIFIER 
INTEGER  DESTOD 

THE  DESTINATION  NODE  OF  AN  OD  PAIR 

INTEGER  ARC  j 

DO  LOOP  INDEX  FOR  ARCS 
INTEGER  PATH 

A  PATH  INDEX 
INTEGER  NUMLIST 

TOTAL  NUMBER  OF  ACTIVE  PATHS  FOR  OD  PAIR  UNDER  CONSIDERATION 
INTEGER  ITER 

SPECIFIC  ITERATION 
INTEGER  N1.N2 

TEMPORARY  VARIABLES 
REAL  MINFDER 

THE  LENGTH  FOR  A  SHORTEST  PATH 
REAL  MINSDER 

THE  SECOND  DERIVATIVE  LENGTH  FOR  THE  SHORTEST  PATH 
REAL  TMINSDER 

TEMPORARY  VALUE  FOR  SECOND  DERIVATIVE  LENGTH  OF  SHORTEST  PATH 
REAL  I NCR 

TOTAL  SHIFT  OF  FLOW  TO  THE  MINIMUM  FIRST  DERIVATIVE  LENGTH  PATH 
REAL  PATH  I NCR  1 

SHIFT  OF  FLOW  FOR  A  GIVEN  PATH  I 
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REAL  FLOW 

C  FLOW  FOR  A  PATH 

REAL  FDER 

C  THE  ACCRUED  LENGTH  ALONG  A  PATH 

REAL  SDER 

C  THE  ACCRUED  SECOND  DERIVATIVE  LENGTH  ALONG  A  PATH 

REAL  TEMP ERROR 

C  TEMPORARY  STORAGE  FOR  CONVERGENCE  ERROR 

REAL  FDLENGTH (NMAXITER) 

C  ARRAY  OF  LENGTHS  OF  PATHS  FOR  AN  OD  PAIR 

REAL  SDLENGTH (NMAXITER) 

C  ARRAY  OF  SECOND  DERIVATIVE  LENGTHS  OF  PATHS  FOR  AN  OD  PAIR 

INTEGER  PATHL I ST (NMAXITER) 

C  ARRAY  OF  ACTIVE  PATHS  FOR  AN  OD  PAIR 

INTEGER  COMMODITY 

C  DO  LOOP  INDEX  FOR  THE  OD  PAIR  ORIGINS 

INTEGER  ORIGIN 

C  SPECIFIC  ORIGIN 

INTEGER  I 

C  DO  LOOP  INDEX 

INTEGER  OD 

C  OD  DO  LOOP  INDEX 

INTEGER  K 

C  DO  LOOP  INDEX 

INTEGER  SHORTEST 

C  THE  SHORTEST  PATH 

LOGICAL  MEMBER  (NNA) 

C  LOGICAL  FOR  AN  ARC  INCLUDED  IN  THE  SHORTEST  PATH 


REAL 

DLENGTH 

c 

REAL 

DIFFERENCE 

D1CAL 

IN 

PATH  LENGTHS  FOR  THE  TRAFFIC 

c 

REAL 

ARC  LENGTH 
D2CAL 

c 

DERIVATIVE 

OF 

ARC  LENGTH 

C 

C  **********************  EXECUTABLE  CODE  ***************************** 

C 

Q  ***************************************** 

C  *  INITIALIZATION 

Q  ***************************************** 

C 

DO  5  ARC=1 ,  NA 
FA (ARC) =0.0 
CONTINUE 


DO  1=1 ,NUMODPAIR 

FP (I) =INPUT_FLOW (I) 

ENDDO 

STARTOD  (NUMCOMMOD+ 1 )  =NUMODPAIR+ 1 

NUMPATH=0 

NUMITER=1 

DO  100  COMMODI TY = 1 , NUMCOMMOD 
OR  I GIN=ORGID (COMMODITY) 

CALL  SP (ORIGIN, COMMODITY) 

DO  10  1=1. NN 

PRED ( I. 1. COMMODITY) =PA(I) 
0  CONTINUE 


LOOP  OVER  OD  PAIRS  OF  COMMODITY 


'  1 


N1=START0D (COMMODITY) 

N2=STARTOD (COMMODITY* 1) -1 
DO  50  OD-N1 ,  N2 

NUMPATH -NUMPATH  1 
PATH ID (NUMPATH) =1 
NEXTPATH (NUMPATH) =0 
FLOW=FP (NUMPATH) 

NODE=DEST (OD) 

DO  WHILE  (NODE. NE. ORIGIN) 
ARC=PA (NODE) 

FA (ARC) =FA (ARC) +FLOW 
NOD  E = S  TAR  TNODE (ARC) 

END  DO 
CONTINUE 
CONTINUE 

INITIALIZE  THE  MEMBER  ARRAY 

DO  70  ARC-1 ,NA 

MEMBER (ARC) = . FALSE . 

CONTINUE 

INITIALIZE  THE  TOTAL  DELAY 
CALL  DELAY (DTOT (NUMITER) ) 

OUTPUT  THE  CURRENT  INFORMATION  TO  DISK 
CALL  PRFLOW 


nonnn  non  noo 
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ARC=PA(NODE) 

IF  (ARC . NE . PRED (NODE , 1 #  COMMODI TY) )  GO  TO  180 
NODE=STARTNODL (ARC) 

END  DO 
GO  TO  500 
END  IF 
C 

180  CONTINUE 

MARK  THE  ARCS  OF  THE  SHORTEST  PATH 

DESTOD=DEST (OD) 

NODE=DESTOD 

DO  WHILE  (NODE. NE. ORIGIN) 

ARC=PA (NODE) 

MEMBER (ARC) = . TRUE . 

NODE=STARTNODE (ARC) 

END  DO 

GENERATE  LIST  OF  ACTIVE  PATHS  FOR  OD  PAIR 


NUMLIST=1 
PATHLI ST ( 1 ) =OD 
PATH=NEXTPATH (OD) 

DO  WHILE  (PATH.GT.O) 

NUML I ST=NUML  I ST+ 1 
PATHLI ST (NUML I ST) =PATH 
PATH=NEXTPATH (PATH) 

END  DO 

DETERMINE  1ST  &  2ND  DERIVATIVE  LENGTH  OF  ACTIVE  PATHS 
ALSO  DETERMINE  WHETHER  THE  CALCULATED  SHORTEST  PATH 
IS  ALREADY  IN  THE  LIST 

SPNEW= . TRUE . 

DO  200  K=1 .NUMLIST 
SAME= . TRUE . 

FDER=0 
SDER=0 
TMINSDER=0 
PATH=PATHLIST (K) 

ITER=PATHID (PATH) 

NODE=DESTOD 

DO  WHILE  (NODE. NE. ORIGIN) 

ARC=PRED (NODE , I TER , COMMODITY) 

CALL  DERIVS (COMMODITY, FA  (ARC) , ARC , D1CAL , D2CAL) 
FDER=FDER+D1CAL 
IF  (.NOT. MEMBER (ARC))  THEN 
SDER=SDER+D2CAL 
SAME=. FALSE. 

ELSE 

SDER=SDER-D2CAL 
TMINSDER=TMINSDER+D2CAL 
END  IF 

NODE =STARIHODE (ARC) 

END  DO 

IF  (SAME)  THEN 
SPNEW=. FALSE. 

SHORTEST=PATH 
_ FDLENGTH (K) =FDER 


-A & 
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MI NFDER-FDER 
MINSDER=TMINSDER 

ELSE 

FDLENGTH(K) =FDER 
SDLENGTH (K) =  SDER 
END  IF 

200  CONTINUE 

***  INSERT  SHORTEST  PATH  IN  PATH  LIST  IF  IT  IS  NEW  *** 

IF  (SPNEW)  THEN 

NUMP ATH=NUMP ATH  + 1 
SHORTEST=NUMPATH 
PATHID (NUMP ATH) =NUMI TER 
NEXTPATH (PATHLIST (NUMLI ST) ) =NUMPATH 
NEXTPATH (NUMP ATH) = 0 
MINFDER=0 
MINSDER=0 
NODE=DESTOD 

DO  WHILE  (NODE. NE. ORIGIN) 

ARC=PA (NODE) 

CALL  DERI VS (COMMODITY. FA (ARC) , ARC . D1CAL , D2CAL) 
MINFDER=MINFDER+D1CAL 
MINSDER=MINSDER+D2CAL 
NODE=STARTNODE (ARC) 

END  DO 
END  IF 

****  UPDATE  PATH  &  LINK  FLOWS  **** 


INCR=0 

TEMPERROR=0 

DO  250  K=l. NUMLI ST 

DLENGTH=FDLENGTH (K) -MINFDER 
IF  (DLENGTH.GT.O)  THEN 
PATH=PATHL I ST (K) 

FLOW=FP (PATH) 

IF  ( (FLOW. EQ. 0.0) .AND. (K.GT.l) )  THEN 
NEXTPATH (PATHLIST (K-l) ) =NEXTP ATH (PATH) 

GO  TO  250 
END  IF 

PATHINCR=DLENGTH/ (SDLENGTH (K) +MINSDER) 

IF  (FLOW. LE.PATHI NCR)  THEN 
FP (PATH) =0 . 0 
PATHINCR=FLOW 
ELSE 

FP (PATH) =FLOW-PATHINCR 
END  IF 

INCR=INCR +PATHI NCR 

TEMPERROR=TEMPERROR+FLOW*DLENGTH/FDLENGTH(K) 
ITER=PATHID (PATH) 

NODE=DESTOD 

DO  WHILE  (NODE. NE. ORIGIN) 

ARC=PRED (NODE , I TER , COMMODI TY) 

FA  (ARC)  =FA  (ARC)  -PATHINCR 
NODE=STARTNODE (ARC) 

END  DO 
END  IF 

250  CONTINUE 


I 


c 

C  ***  UPDATE  THE  ERROR  CRITERION  *** 

C 

CURERROR=AMAXl ( CUR  ERROR , TEMP  ERROR/ 1 NPUT_F LOW (OD) ) 
C 

c  ****  UPDATE  FLOWS  FOR  SHORTEST  PATH  **** 

C 

FP (SHORTEST) =FP (SHORTEST) + INCR 
NODE=DESTOD 

DO  WHILE  (NODE. NE. ORIGIN) 

ARC~PA (NODE) 

FA (ARC) =FA (ARC) ♦ INCR 
MEMBER (ARC) = . FALSE . 

NODE = S  TAR  TNODE (ARC) 

END  DO 
C 

500  CONTINUE 

C 

c  *****  END  OF  LOOP  FOR  OD  PAIRS  CORRESPONDING  TO  COMMODITY 

C  *****  UPDATE  TOTAL  DELAY 

C 

CALL  DELAY (DTOT(NUMI TER) ) 

C 

1000  CONTINUE 

C 

C  CHECK  IF  THE  #  OF  ACTIVE  PATHS  EXCEED  THE  ALLOCATED  NUMBER 

C 

IF  (NUMPATH .  GT .  NNUMPATH)  THEN 

WRITE  (6, *) 'MAX  #  OF  ALLOCATED  PATHS  EXCEEDED’ 

STOP 
END  IF 
C 

C  OUTPUT  THE  CURRENT  SOLUTION  TO  DISK 

C 

CALL  PRFLOW 
C 

C  *****  END  OF  ITERATION  ***** 

C 

C  ***  IF  THE  ERROR  IS  SMALLER  THAN  TOL,  OR  THE  LIMIT  ON 

C  THE  NUMBER  OF  ITERATIONS  IS  REACHED  RETURN 

C  ELSE  GO  FOR  ANOTHER  ITERATION 

C 

IF  ( (CURERROR.LT.  TOL)  .OR.  (NUMI  TER.  EQ.  MAXI  TER))  THEN 
RETURN 

ELSE 

GO  TO  110 
END  IF 
C 

END 

c  *.****«*.*.**.  END  OF  MULTIFLO  **************** 
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SHORTHEAP 

’  SHORTHEAP '  SOLVES  THE  SHORTEST  PATH  PROBLEM  BY 
DIJKSTRA'S  ALGORITHM  AND  A  HEAP  DATA  STRUCTURE. 

THIS  ALGORITHM  SHOULD  BE  USED  WHEN  THE  NUMBER  OF 
DESTINATIONS  FOR  EACH  COMMODITY  IS  SMALL  RELATIVE 
TO  THE  TOTAL  NUMBER  OF  NODES. 

INPUT: 

S  -  THE  STARTING  NODE 

COMMODITY  -  THE  CORRESPONDING  COMMODITY 
OUTPUT: 

PA  (I)  -  THE  LAST  ARC  ON  THE  SHORTEST  PATH  ENDINC  AT  NODE  I 
DIST(I)  -  THE  SHORTEST  DISTANCE  TO  NODE  I 


SUBROUTI NE  SP  ( S , COMMODI TY) 

IMPLICIT  NONE 

******************  INCLUDE  COMMON  BLOCKS  ** 

INCLUDE  'PARAM.DIM' 

INCLUDE  'NETWRK.PRM' 

INCLUDE  ’ PATHS. BLK' 

*********.****.**  local  VARIABLE  DEFINITIONS 


REAL 


REAL 


REAL 


INTEGER 


INTEGER 


INTEGER 


INTEGER 


INTEGER 


INTEGER 


INTEGER 


INTEGER 


INTEGER 


INTEGER 


INTEGER 


INTEGER 


MIN 

TEMPORARY  MINIMUM  VALUE 
Dl , D2 , DP 
NODE  DISTANCE 
XLARGE 

BIG  X  BY  DEFAULT 
S 

INPUT  NODE 
COMMODITY 
INPUT  COMMODITY 
P 

NODE  ALONG  THE  PATH  OF  S  TO  DESTINATIONS 
I 

DO  LOOP  INDEX 
J 

DO  LOOP  INDEX 
ARC 

DO  LOOP  INDEX 
ND 

A  NODE  INDEX 
DNUMBER 

#  OF  DESTINATIONS  FOR  COMMODITY 
N1 

TEMPORARY  VARIABLE 
N2 

TEMPORARY  VARIABLE 
UPNODE , DOWNNODE , DOWNNODE1 , LASTNODE 
VARIABLES  USED  IN  UPDATING  THE  HEAP  ARRAY 
CURRANK , NEWRANK 
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C  VARIABLES  USED  IN  UPDATING  THE  HEAP  ARRAY 

INTEGER  END: ZAP 

C  MARKS  THE  LAST  ELEMENT  OF  THE  HEAP  ARRAY 

INTEGER  RANK (NNN) 

C  RANK (NODE)  GIVES  THE  RANK  OF  NODE  IN  THE  HEAP 

INTEGER  NRANK(NNN) 

C  NRANX(I)  GIVES  THE  NODE  OF  RANK  I  IN  THE  HEAP 

REAL  D1CAL 

C  FIRST  DERIVATIVE  OF  DELAY  WITH  RESPECT  TO  LOAD 

LOGICAL  FIRSTITER 

C  TRUE  IF  THIS  IS  THE  FIRST  ITERATION 

LOGICAL  SCAN (NNN) 

C  LOGICAL  INDICATING  THAT  A  NODE  HAS  BEEN  SCANNED 

LOGICAL  DSTATUS (NNN) 

LOGICAL  SPECIFYING  IF  A  NODE  IS  A  DESTINATION 

******************  EXECUTABLE  CODE  ********************* 

XLARGE=1E1S 
D1CAL=1 .0 
P=S 

DO  10  1=1, NN 

DI ST ( I ) =XLARGE 
SCAN (I ) = .FALSE . 

DSTATUS (I) =. FALSE. 

D  CONTINUE 

DIST (S) =0 

IF  (NUMITER . EQ. 1)  THEN 
FIRSTITER=.TRUE. 

ELSE 

FIRSTITER=. FALSE. 

END  IF 

MARK  THE  DESTINATION  NODES 

Nl=STARTOD (COMMODITY) 

N2= STAR TOD (COMMODI TY  + 1 ) -1 
DNUMBER-N2-N1+ 1 
DO  15  I=N1.N2 

DSTATUS (DEST ( I ) ) = . TRUE . 

5  CONTINUE 

INITIALIZE  THE  HEAP  FLOOR 

ENDHEAP=0 

*****  SCAN  NODE  P  ***** 

1000  CONTINUE 

SCAN (P) =.TRUE . 

IF  (DSTATUS (P))  THEN 

IF  (DNUMBER . EQ . 1 )  RETURN 
DNUMBER=DNUMBER- 1 
END  IF 

IF  (FRSTOU(P) .NE.O)  THEN 
DP=DIST(P) 

DO  20  ARC=FRSTOU (P) , LASTOU (P) 

ND=ENDNODE (ARC) 

IF  (. NOT. SCAN (ND))  THEN 

IF  (.NOT. FIRSTITER)  THEN 


*«***-«* 
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CALL  DERI VI (COMMODITY. FA (ARC) , ARC.D1CAL) 
END  IF 
D2=DIST (ND) 

C  IF  ND  HAS  NOT  BEEN  LABELLED  INSERT  IT  IN  THE  HEAP 

IF  (D2.EQ.XLARGE)  THEN 
ENDHEAP = ENDHEAP  + 1 
RANK (ND) = ENDHEAP 
NRANK (ENDHEAP) =ND 
END  IF 
D1=DP+D1CAL 
IF  (D1.LT.D2)  THEN 
PA (ND) =ARC 
DIST (ND) =D1 
CURRANK=RANK (ND) 

SO  NEWRANK=INT (CURRANK/2) 

IF  (NEWRANK . GE . 1 )  THEN 
UPNODE=NRANK (NEWRANK) 

IF  (D1 . LT .DIST (UPNODE) )  THEN 
NRANK (CURRANK) =UPNODE 
RANK (UPNODE) =CURRANK 
CURRANK=NEWRANK 
GO  TO  50 
END  IF 
END  IF 

NRANK (CURRANK) =ND 
RANK (ND) =CURRANK 
END  IF 
END  IF 

20  CONTINUE 

END  IF 

*******  FIND  next  NODE  TO  SCAN  ******* 

TEST  FOR  ERROR 
IF  (ENDHEAP.EQ.O)  THEN 

WRITE (6,*)  'ERROR  IN  THE  SHORTEST  PATH  POUTINE' 
STOP 
END  IF 
P=NRANK  (1) 

RESTRUCTURE  HEAP  ARRAYS 

LASTNODE=NRANK  (ENDHEAP) 

ENDHEAP = ENDHEAP - 1 
D1=DI ST (LASTNODE) 

CURRANK=1 
100  NEWRANK=CURRANK + CURRANK 

IF  (NEWRANK. LE. ENDHEAP)  THEN 
DOWNNODE=NRANK (NEWRANK) 

IF  (NEWRANK. EQ. ENDHEAP)  THEN 
DOWNNODE l=DOWNNODE 
ELSE 

DOWNNODE 1 =NRANK (NEWRANK+1) 

END  IF 

IF  (DIST (DOWNNODE) .LE. DIST (DOWNNODE1) )  THEN 
IF  (Dl.GT. DIST (DOWNNODE))  THEN 
NRANK (CURRANK) = DOWNNODE 
RANK  (DOWNNODE)  =  CURRANK 
CURRANK=NEWRANK 
GO  TO  100 


END  IF 
ELSE 

IF  (D1 .GT.DIST (DOWNNODE1) )  THEN 
NRANK (CURRANK) =DOWNNODEl 
RANK (DOWNNODE1) = CURRANK 
CURRANK=NEWRANK+ 1 
GO  TO  100 
END  IF 
END  IF 
END  IF 

NRANK (CURRANK) =LASTNODE 
RANK (LASTNODE) = CURRANK 
GO  TO  1000 

END 
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c 

C  SHORTPAPE 

C  'SHORTPAPE'  SOLVES  THE  SHORTEST  PATH  PROBLEM  BY 

C  PAPE’S  MODIFICATION  OF  BELLMAN'S  ALGORITHM. 

C 

C  INPUT: 

C  S  -  THE  STARTING  NODE 

C  COMMODITY  -  THE  CORRESPONDING  COMMODITY 

C 

C  OUTPUT: 

C  PA (I)  -  THE  LAST  ARC  ON  THE  SHORTEST  PATH  ENDING  AT  NODE  I 

C  DIST(I)  -  THE  SHORTEST  DISTANCE  TO  NODE  I 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

SUBROUTINE  SP (S , COMMODITY) 

C 


IMPLICIT  NONE 


******************  INCLUDE  COMMON  BLOCKS 


INCLUDE  'PARAM.DIM' 

INCLUDE  'NETWRK.PRM' 

INCLUDE  ' PATHS. BLK’ 

*****************  LOCAL  VARIABLE  DEFINITIONS  ******* 

REAL  Dl.DP 

NODE  DISTANCE 
REAL  X LARGE 

BIG  X  BY  DEFAULT 
INTEGER  I LARGE 

INTEGER  LARGER  THAN  THE  NUMBER  OF  NODES 
INTEGER  S 

INPUT  NODE 
INTEGER  COMMODITY 

INPUT  COMMODITY 
INTEGER  P 

NODE  PRESENTLY  SCANNED 
INTEGER  I 

DO  LOOP  INDEX 
INTEGER  ARC 

EX)  LOOP  INDEX 
INTEGER  ND 

A  NODE  INDEX 
INTEGER  N1 

TEMPORARY  VARIABLE 
INTEGER  N2 

TEMPORARY  VARIABLE 
INTEGER  ENDQUEUE 

MARKS  THE  LAST  ELEMENT  OF  THE  QUEUE  ARRAY 
REAL  D1CAL 

FIRST  DERIVATIVE  OF  DELAY  WITH  RESPECT  TO  FLOW 
LOGICAL  FIRSTITER 

TRUE  IF  THIS  IS  THE  FIRST  ITERATION 
INTEGER  Q(NNN) 

QUEUE  OF  NODES  TO  BE  SCANNED 

*********.*.•*•«•*  EXECUTABLE  CODE  **************** 


c 


XLARCE=1E15 
ILARGE^NNN+l 
D1CAL=1 .0 
DO  10  1=1, NN 

DIST (I) =XLARGE 

Q(i)=o 

10  CONTINUE 

IF  (NUMITER.EQ.l)  THEN 
F IRSTI TER= . TRUE . 

ELSE 

F IRSTI TER=. FALSE. 

END  IF 
DIST (S) =0 
Q(S)=ILARGE 
ENDQUEUE=S 
P=S 
C 

C  *********  START  OF  MAIN  ALGORITHM  ********* 

C 

100  CONTINUE 

C 

C  *****  SCAN  NODE  P  ***** 

C 

Nl=FRSTOU(P) 

IF  (N1.EQ.0)  GO  TO  201 
N2=LASTOU (P) 

DP=DI ST (P) 

DO  200  ARC=N1,N2 
ND=ENDNODE (ARC) 

IF  (. NOT. F IRSTI TER)  THEN 

CALL  DERIV1 (COMMODITY, FA{ARC) , ARC, D1CAL) 

END  IF 
D1=DP+D1CAL 

C  »“  IF  NO  IMPROVEMENT’  TAKE  ANOTHER  ARC  *** 

IF  (Dl.GE.DIST(ND))  GO  TO  200 
C  ***  CHANGE  DISTANCE  AND  LABEL  OF  NODE  ND  *  *  * 

PA  (ID)  “ARC 
DIST (ND) =D1 

IF  (Q (ND) )  160,140,200 

C  *»«  IF  ND  HAS  NEVER  BEEN  SCANNED  INSERT  IT  AT  THE  END 

C  OF  THE  QUEUE  *** 

140  Q (ENDQUEUE) =ND 

ENDQUEUE=ND 
Q  (ND) = I LARGE 
GO  TO  200 

C  ***  IF  ND  HAS  ALREADY  BEEN  SCANNED  ADD  IT  AT  THE 

C  BEGINNING  OF  THE  QUEUE  AFTER  NODE  P  *** 

160  Q(ND)=Q(P) 

Q(P)  =ND 

IF  (ENDQUEUE . EQ.P)  ENDQUEUE=ND 

200  CONTINUE 
C 

C  ***  GET  NEXT  NODE  FROM  THE  TOP  OF  THE  QUEUE  *** 

C 

201  Nl=Q (P) 

C 

C  ***  FLAG  P  AS  HAVING  BEEN  SCANNED  *** 

C 


Q(P)=-1 


P-Nl 


***  IF  THE  QUEUE  IS  NOT  EMPTY  GO  BACK  TO  SCAN  NEXT  NODE 

IF  (P.LT.ILARGE)  GO  TO  100 

RETURN 

END 
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c 

C  DELAY 

c 

C  DELAY  COMPUTES  THE  TOTAL  M/M/1  DELAY  IN  ROUTING  COMMODITIES  FROM 

C  SOURCES  TO  SINKS. 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

SUBROUTINE  DELAY (DT) 

IMPLICIT  NONE 

*****************  INCLUDE  COMMON  BLOCKS  ************************** 

INCLUDE  'PARAM.DIM' 

INCLUDE  ’ PATHS. BLK’ 

INCLUDE  * NETWRK.PRM* 

INCLUDE  ' CONVRG . PRM * 

********************  ARGUMENT  DEFINITIONS  ************************* 

ON  OUTPUT: 

REAL  DT 

TOTAL  SYSTEM  DELAY 

******************  EXTERNAL  FUNCTIONS  REFERENCED  ****************** 

REAL  DCAL 

DELAY  AS  A  FUNCTION  OF  FLOW 

******************  LOCAL  VARIABLE  DEFINITIONS  ********************** 

INTEGER  K 

DO  LOOP  INDEX 

***********************  EXECUTABLE  CODE  **************************** 

LOOP  OVER  ALL  LINKS  AND  ACCRUE  TOTAL  DELAY 
DT=0 . 

DO  50  K=1,NA  _ 

DT=DT+DCAL (FA (K) , K) 

50  CONTINUE 
C 


RETURN 

END 
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DCAL 

'DCAL'  COMPUTES  THE  DELAY  ACROSS  A  SPECIFIED  ARC  GIVEN  THE  FLOW. 

THE  DELAY  IS  ASSUMED  TO  BE  CONSISTENT  WITH  M/M/1  QUEUEING  FOR 
FLOWS  BELOW  A  MAXIMUM  UTILIZATION  AND  QUADRATIC  BEYOND  WITH 
CONTINUITY  IN  THE  DERIVATIVES  AT  THE  MAXIMUM  UTILIZATION. 

:ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

REAL  FUNCTION  DCAL (X, ARC) 

IMPLICIT  NONE 


INCLUDE  COMMON  BLOCKS 


INCLUDE  ’PARAM.DIM' 

INCLUDE  'NETWRK.PRM' 

INCLUDE  ' CONVRG . PRM ' 

INCLUDE  ' PATHS . BLK ' 

*****************  ARGUMENT  DEFINITIONS 


REAL  X 

INPUT  FLOW  FOR  THE  ARC 
INTEGER  ARC 

INPUT  ARC 

****************  LOCAL  VARIABLE  DEFINITIONS  ************** 

REAL  RATE 

MAXIMUM  LINK  CAPACITY 
REAL  Y 

TEMPORARY  VARIABLE 
REAL  Z 

TEMPORARY  VARIABLE 
REAL  QO 

ZEROTH  ORDER  TERM  IN  THE  QUADRATIC  APPROXIMATION  FOR 
OVERLOADED  LINKS 
REAL  Q1 

FIRST  ORDER  TERM  IN  THE  QUADRATIC  APPROXIMATION 
REAL  Q2 

SECOND  ORDER  TERM  IN  THE  QUADRATIC  APPROXIMATION 
REAL  EXCESS 

FLOW  BEYOND  THE  MAXIMUM  ALLOWABLE  UTILIZATION 

**********************  EXECUTABLE  CODE  ****** ** *••*••*••* * 

RATE=BITRATE  (ARC) 

Y=MAXUTI ‘RATE 

M/M/1  DELAY 

IF(X.LT.Y)  THEN 

DCAL=X/ (RATE-X) 

ELSE 


QUADRATIC  APPROXIMATION  TO  AVOID  OVERFLOWS 


EXCESS=X-Y 


nonn  o  nnnn  n  nn  nnnn  nno  nnnnnonnoon 


Z=RATE-Y 

QO=Y/Z 

Q1=Q0/ (MAXUTI *Z) 

Q2=Q1/Z 

DCAL=Q0+Q1*EXCESS+Q2* EXCESS*  *2 
ENDIF 
RETURN 
END 

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DERI  VS 

'DERI VS'  COMPUTES  THE  DERIVATIVES  OF  DELAY  WITH  RESPECT  TO  FLOW  FOR 
LINKS.  BELOW  A  MAXIMUM  UTILIZATION.  M/M/1  DELAY  IS  ASSUMED  TO  APPLY 
WHEREAS  A  QUADRATIC  APPROXIMATION  IS  ASSUMED  FOR  UTILIZATIONS  BEYOND 
THE  MAXIMUM.  THE  DERIVATIVES  ARE  CONTINUOUS  AT  THE  MAXIMUM 
UTILIZATION. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTI NE  DER I VS (COMMOD I TY . X . ARC , D 1 CAL , D2CAL) 

IMPLICIT  NONE 

********************  INCLUDE  COMMON  BLOCKS  **•••**••*•***•*****•*** 

INCLUDE  ' PARAM.DIM' 

INCLUDE  'NETWRK.PRM' 

INCLUDE  ' CONVRG . PRM ' 

INCLUDE  ' PATHS . BLK ' 

*******************  ARGUMENT  DEFINITIONS  ************************** 

ON  INPUT: 

INTEGER  COMMODITY 

THE  CORRESPONDING  COMMODITY 

REAL  X 

FLOW  IN  THE  SPECIFIED  LINK 
INTEGER  ARC 

THE  SPECIFIED  LINK 


ON  OUTPUT: 

REAL  D1CAL 

ARC  LENCTH  (1ST  DERIVATIVE  OF  DELAY) 
REAL  D2CAL 

FIRST  DERIVATIVE  OF  ARC  LENGTH 
****************  LOCAL  VARIABLE  DEFINITIONS 


C 

C 

c 

c 

c 


REAL  ■  MAXI 

MAXIMUM  ALLOWABLE  FLOW  FOR  LINK  FOR  M/M/1  QUEUEING  DELAY 
REAL  RATE 

THE  MAXIMUM  FLOW  CAPACITY  FOR  THE  LINK 
REAL  EXCESS 

FLOW  BEYOND  THE  MAXIMUM  ALLOWABLE  FLOW 
REAL  D1 

TEMPORARY  VARIABLE 
REAL  T 

TEMPORARY  VARIABLE 


non  n  nn  nnnnn  non  onnnnnnnnoon  non  nnn  o  nnn 


T 


********************  EXECUTABLE  code 


RATE=BI TRATE (ARC) 

MAXI=MAXUTI ‘RATE 
EXCESS=X-MAXI 

IF  (EXCESS. LE. 0.0)  THEN 

DERIVATIVES  OF  M/M/1  QUEUEING  DELAY 


T=RATE-X 
D1CAL=RATE/T*  *  2 
D2CAL=2 . 0  *D1CAL/T 

ELSE 


DERIVATIVES  OF  THE  QUADRATIC  APPROXIMATION 

T=RATE-MAXI 
D1=RATE/T**2 
D2CAL=2.0*D1/T 
D1CAL=D1+D2CAL* EXCESS 
END  IF 
RETURN 
END 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

DERIV1 

'DERIV1 *  COMPUTES  THE  FIRST  DERIVATIVE  OF  DELAY  WITH  RESPECT 
TO  FLOW  FOR  LINKS.  BELOW  A  MAXIMUM  UTILIZATION,  M/M/1  DELAY  IS 
ASSUMED  TO  APPLY  WHEREAS  A  QUADRATIC  APPROXIMATION  IS  ASSUMED  FOR 
UTILIZATIONS  BEYOND  THE  MAXIMUM.  THE  DERIVATIVES  ARE  CONTINUOUS 
AT  THE  MAXIMUM  UTILIZATION. 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

SUBROUTINE  DERIV1 (COMMODITY, X, ARC, D1CAL) 

IMPLICIT  NONE 

********************  INCLUDE  COMMON  BLOCKS  a*********************** 


INCLUDE  'PARAM.DIM' 

INCLUDE  ' NETWRK . PRM ' 

INCLUDE  ' CONVRG . PRM ' 

INCLUDE  ' PATHS. BLK’ 

*******************  ARGUMENT  DEFINITIONS 


ON  INPUT: 

INTEGER  COMMODITY 

THE  CORRESPONDING  COMMODITY 

REAL  X 

FLOW  IN  THE  SPECIFIED  LINK 
INTEGER  ARC 

THE  SPECIFIED  ARC 

ON  OUTPUT: 


k 


REAL  D1CAL 

ARC  LENGTH  (1ST  DERIVATIVE  OK  DELAY) 

.  local  variable  definitions 


REAL  MAXI 

MAXIMUM  ALLOWABLE  FLOW  FOR  LINK  FOR  M/M/1  QUEUEING  DELAY 
REAL  RATE 

THE  MAXIMUM  FLOW  CAPACITY  FOR  THE  LINK 
REAL  EXCESS 

FLOW  BEYOND  THE  MAXIMUM  ALLOWABLE  FLOW 
REAL  D1 

TEMPOILARY  VARIABLE 
REAL  T 

TEMPORARY  VARIABLE 
REAL  D2CAL 

TEMPORARY  VARIABLE 


RATE- B I  IRATE (A 
MAXI  MAXUIV  »?A 

EXCESS  x-  max: 


IVRaTE  X 
D^CAD-RATt'/I  **2 


DLR XVA'i  IVE  Of  'II «  QUADRAT!  l  APP1<’ •  >X  MA *  1  ^ ■ 


nnn  nnnn  n  non  non  nnnonnoon 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

PRFLOW 

' PRFLOW ’  OUTPUTS  INTERMEDIATE  RESULTS  IN  THE  MULTIFLO  ALGORITHM. 
ITERATION  #,  DELAY,  NUMBER  OF  ACTIVE  PATHS  GENERATED  AND 
CONVERGENCE  ARE  THE  PRIMARY  OUTPUTS. 

xccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  PRFLOW 
IMPLICIT  NONE 

****************  INCLUDE  COMMON  BLOCKS  *******4** ***************** 

INCLUDE  'PARAM.DIM* 

INCLUDE  'NETWRK.PRM'  - 

I NCLUDE  ' CONVRG . PRM ' 

INCLUDE  ' PATHS . BLK ' 

**************  LOCAL  VARIABLE  DEFINITIONS  *********************** 

LOGICAL  FIRFLG 

FIRST  PASS  FLAG  FOR  OUTPUT  CONTROL 
INTEGER  I 

DO  LOOP  INDEX 

*****************  LOCAL  DATA  INITIALIZATION  ******************* 


DATA  FIRFLG/. TRUE./ 


ON  THE  VERY  FIRST  PASS,  OUTPUT  THE  CONTENTS  OF  INPUT  BLOCKS  TO  FILE 


IF  (FIRFLG)  THEN 


WRITE (6, 
WRITE (6, 


WRITE (6, 

WRITE (6, 

WRITE (6,*) '**«** ************************************ 

INITIALIZATION  DATA 


WRITE (6, 
WRITE (6, 


***************************************** 

*  MULTIFLO  SUMMARY  *' 

A**************************************** 


***************************************** 


) 'NETWORK  SPECIFICATION  DATA:’ 
)’  ' 

) 'NODE  SPECIFICATIONS' 

) 'NUMBER  OF  NODES:', NN 
)  'NODE  #  FRSTOU 


WRITE (6, 

WRITE (6 , 

WRITE (6, 

WRITE (6 , 

WRITE (6, 

WRITE (6, 

DO  1=1, NN 

WRITE (6 , *) I , FRSTOU (I ) ,LASTOU(I) 
END  DO 

WRITE (6,*)'  ' 

WRITE (6,*) 'LINK  SPECIFICATIONS:' 
WRITE (6,*) 'NUMBER  OF  LINKS: ' ,NA 
WRITE (6, *) 'LINK  #  STARTNODE 

DO  1=1, NA 


LASTOU' 


ENDNODE 


WRITE (6,*) I, STARTNODE (I) , ENDNODE (I) , BITRATE (I) 
END  DO 

WRITE (6,*)’  ' 

WRITE (6.*) 'COMMODITY  SPECIFICATIONS' 

WRITE (6,*) 'NUMBER  OF  COMMODITIES :* ,NUMCOMMOD 


BITRATE ' 


WRITE (6, ‘)  ' COMMOD  #  ORG I D  STARTOD * 

DO  1=1 .NUMCOMMOD 

WRITE (6 ,  * ) I ,ORGID (I ) ,STARTOD(I) 

END  DO 

WRITE (6 ,  * ) *  ' 

WRITE (6, *) 'OD  PAIR  SPECIFICATIONS'  — 

WRITE  (6,*) 'NUMBER  OF  OD  PAIRS:  \NUMODPAIR 

WRITE  (6 ,  *)  'OD  PAIR  #  DEST  INPUT  FLOW' 

DO  1=1, NUMODPAIR 

WRITE (6. *) I ,DEST (I) , INPUT_FLOW (I) 

END  DO 

WRITE (6,*)'  ' 

WRITE (6, *)  ****************************************** 
WRITE (6,*)'*  MULTIFLO  DATA  BY  ITERATION  * ' 

WRITE (6, *)’****************************•***********' 
WRITE (6,*) ’ITERATION  #  TOTAL  DELAY  CONVERGENCE 

WRITE (6,*)'  ERROR 

WRITE (6,  *)  ’ 

FIRFLG=. FALSE. 

END  IF 

IF  (NUMITER.GT.O)  THEN 

WRI TE (6 , * ) NUMI TER . DTOT (NUMI TER) , CUR ERROR , NUMPATH 
END  IF 
RETURN 
END 


NUMBER  OF' 
ACTIVE ' 
PATHS' 


nnn  n  n  o  o  o  n  nnonno 


’INCLUDE’  FILE  P ARAM. DIM 


’PARAM.DIM' 

************ 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 

PARAMETER 


CONTAINS  THE  ARRAY  DIMENSIONS 
*******  NETWORK  PARAMETERS  * 


NNN=100 

MAXIMUM  NUMBER  OF  NODES 
NNA=500 

MAXIMUM  NUMBER  OF  ARCS 
NNUMOD=1000 

MAXIMUM  NUMBER  OF  OD  PAIRS 
NNUMPATH=10000 

MAXIMUM  NUMBER  OF  PATHS  FOR  CONSIDERATION 
NMAXITER=50 

MAXIMUM  NUMBER  OF  ITERATIONS  ALLOWED 
NNORIG=100 

MAXIMUM  NUMBER  OF  COMMODITIES 
NINDEX=100000 

MAXIMUM  NUMBER  OF  ELEMENTS  OF  PATH 
DESCRIPTION  ARRAY  (USED  IN  MULTIFLOl) 


C->  C"* 


1 


u./ 


c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


' INCLUDE  *  FILE  NETWRK.PRM 

’NETWRK.PRM'  CONTAINS  THE  NETWORK  SPECIFICATION  PARAMETERS 

COMMON  /NETWORK/ 

NN .  FRSTOU , LASTOU , 

NA,  STARTNODE , ENDNODE , BITRATE , 

NUMCOMMOD , ORGID , STARTOD , 

NUMODPAIR ,DEST, INPUT_FLOW 


INTEGER*  2 
INTEGER* 2 
INTEGER *2 


NN 

NUMBER  OF  NODES  IN  THE  NETWORK 
FRSTOU (NNN) 

THE  FIRST  ARC  EMANATING  FROM  A  NODE 
LASTOU  (NNN) 

THE  FINAL  ARC  EMANATING  FROM  A  NODE 


INTEGER *2 
INTEGER* 2 
INTEGER *2 
REAL 


NA 

NUMBER  OF  LINKS  (ARCS)  IN  THE  NETWORK 
STARTNODE (NNA) 

THE  START  NODE  FOR  AN  ARC 
ENDNODE (NNA) 

THE  END  NODE  FOR  AN  ARC 
BITRATE (NNA) 

THE  LINK  CAPACITY  IN  BITS/SECOND 


INTEGER* 2 
INTEGER* 2 
INTEGER *2 


NUMCOMMOD 

THE  NUMBER  OF  COMMODITIES  IN  THE  NETWORK 
ORGID (NNORIG) 

THE  NODE  NUMBER  OF  THE  ORIGIN 
STARTOD (NNORIG) 

THE  POINTER  TO  THE  STARTING  NODE  IN  AN  OD  PAIR 


INTEGER* 2 
INTEGER* 2 
REAL 


NUMODPAIR 

THE  NUMBER  OF  OD  PAIRS 
DEST (NNUMOD) 

THE  DESTINATION  NODE  OF  TRAFFIC  IN  AN  OD  PAIR 
INPUT_FLOW (NNUMOD) 

THE  INPUT  TRAFFIC  TO  THE  NODE  IN  BITS/SECOND 


o  n  n  n  n 


'INCLUDE'  FILE  CONVRG.PRM 

' CONVRG . PRM ’  CONTAINS  THE  CONVERGENCE  PARAMETERS  FOR  THE 
NETWORK  FLOW  PROBLEM 

COMMON  /CONVRG/ 

&  MAXITER,TOL,MAXUTI , OUTPFL 

C 

INTEGER  MAX I TER 

C  MAXIMUM  NUMBER  OF  ITERATIONS  IN  THE  SOLUTION 

REAL  TOL 

C  TOLERANCE  ON  SOLUTION  ACCURACY 

REAL  MAXUTI 

C  MAXIMUM  UTILIZATION  FOR  M/M/1  QUEUE  DELAY 

LOGICAL  OUTPFL 

C  OUTPUT  CONTROL  VARIABLE 


i 


C  'INCLUDE*  FILE  PATHS . BLK 

C 

C  'PATHS. ELK'  DEFINES  THE  ARRAYS  NECESSARY  TO  MAINTAIN 

C  PATH  FLOWS  AND  DESCRIPTION. 

C 

COMMON  /PATHS/ 

PA . FA , PATH I D , NEXTPATH , FP , DI ST , DTOT , CURERROR , 
NUMP ATH , NUM ITER 


INTEGER* 2 

PA(NNN) 

c 

THE  LAST  ARC  ON  A  SHORTEST  PATH  TO  A  NODE 

REAL 

FA(NNA) 

c 

THE  FLOW  IN  ANY  GIVEN  LINK  (ARC) 

INTEGER 

PATH I D (NNUMPATH) 

c 

THE  PATH  IDENTIFIER  FOR  ANY  GIVEN  PATH 

INTEGER 

NEXTPATH (NNUMPATH) 

c 

THE  NEXT  PATH  FOR  THE  SAME  OD  PAIR 

REAL 

FP (NNUMPATH) 

c 

THE  FLOW  OF  A  PATH 

REAL 

DIST(NNN) 

c 

SHORTEST  DISTANCE  TO  A  NODE  FROM  THE  ORIGIN 

REAL 

DTOT (NMAXITER) 

c 

THE  TOTAL  DELAY  BY  ITERATION 

INTEGER 

NUMITER 

c 

CURRENT  ITERATION  NUMBER 

REAL 

CURERROR 

c 

CONVERGENCE  ERROR  (NORMALISED  %  OF  FLOW  NOT  ON 

c 

A  SHORTEST  PATH) 

INTEGER 

NUMPATH 

c 

NUMBER  OF  GENERATED  PATHS 

3  <.i  5 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  SETUP 

c 

C  'SETUP’  ACCEPTS  INPUTS  FROM  THE  TERMINAL  AND  CREATES  DATA  SETS 

C  THAT  REPRESENTS  NETWORKS  AND  LOADS  IN  A  FORM  SUITABLE  FOR 

C  PROGRAM  'MULTI FLO' 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

PROGRAM  SETUP 
IMPLICIT  NONE 


C 

c  ********************  INCLUDE  COMMON  BLOCKS  *********************** 

C 

INCLUDE  ' PARAM.DIM' 

INCLUDE  ' NETWRK .PRM ' 

C 

C  ***************  LOCAL  VARIABLE  DEFINITIONS  *********************** 

C 


C 

C 

C 

C 

C 

C 

C 

C 

C 

c 

c 

c 

c 


200 


100 

c 

c 

c 


INTEGER  TERM I NAL_NODE 

THE  END  NODE  OF  A  LINK 
INTEGER  DESTOD 

THE  DESTINATION  NODE  OF  AN  OD  PAIR 
REAL  BPS 

MAXIMUM  LINK  CAPACITY 
INTEGER  NUMARC 

NUMBER  OF  OUTGOING  ARCS  FOR  A  NODE  IN  THE  NETWORK 
REAL  TRAFFIC 

SPECIFIED  INPUT  TO  AN  OD  PAIR 
INTEGER  I 

DO  LOOP  INDEX 
INTEGER  J 

DO  LOOP  INDEX 
INTEGER  NOD 

NUMBER  OF  OD  PAIRS  ASSOCIATED  WITH  A  COMMODITY 

******************  EXECUTABLE  CODE  ****************************** 

GET  THE  NODE  SPECIFICATIONS 


NA=0 

WRITE  (6,*) 'INPUT  THE  #  OF  NODES' 

READ (5, *) NN 
DO  1=1, NN 

WRITE (6,*) ’FOR  NODE*. I,’  ENTER  #  OF  ARCS  EXITING  THE  NODE’ 

READ (5 , * . ERR= 200 ) NUMARC 
IF (NUMARC. GE.O)  THEN 
DO  J=l. NUMARC 

WRITE (6.*) 'FOR  ARC’,J,’  AT  NODE’, I,’  ENTER  TERMINAL  NODE’, 
'  AND  MAXIMUM  BITS/S’ 

ASK  THE  SAME  QUESTION  ON  ERRORS 

READ (5, * , ERR=100) TERMINAL_NODE, BPS 
IF (TERMINAL_NODE .GT.NN)  THEN 

WRITE (6,*) ’TERMINAL  NODE  OUT  OF  BOUNDS’ 

GO  TO  100 


r 


ELSE 


C  ENTER  LINK  BEGIN  AND  END  NODES 

C 

NA-NA+ 1 

ENDNODE (NA) =TERMINAL_NODE 
BITRATE (NA) =  BPS 
END  IF 

STARTNODE (NA) =  1 
END  DO 

FRSTOU ( I ) =NA-NUMARC+ 1 
LASTOU ( I ) =NA 

ELSE 

WRITE (6, *) 'NEGATIVE  ARCS  ILLEGAL' 

GO  TO  200 
END  IF 
END  DO 
C 

C  OD  PAIRS  SETUP 

C 

1000  WRITE (6, *) 'ENTER  THE  NUMBER  OF  COMMODITIES  IN  THE  NETWORK' 

READ (5, * ,ERR=1000)NUMCOMMOD 
NUMODP AIR=0 
DO  1  =  1, NUMCOMMOD 

300  WRITE (6,*) 'ENTER  THE  ORIGIN  ID  AND  NUMBER  OF  DESTINATIONS  FOR  ' 

&  'COMMODITY' , I 

READ (5 ,  *  , ERR=300) ORGID (I ) ,NOD 
IF (ORGID(I) .LE.NN)  THEN 
DO  J-l.NOD 

400  WRITE (6,*) 'ENTER  THE  DESTINATION' , J, '  AND  TRAFFIC  FOR  '. 

&  '  COMMODITY' 

C 

C  ASK  THE  SAME  QUESTION  ON  ERRORS 

C 

READ (5, * ,ERR=400)DESTOD, TRAFFIC 
I F (DESTOD . GT . NN)  THEN 

WRITE  (6, *)  'DESTINATION  OD  OUT  OF  BOUNDS ,  MAXIMUM  '  , NN 
GO  TO  400 

ELSE 

NUMODP A I R=NUMODPA I R  + 1 
DEST (NUMODPAIR) =DESTOD 
INPUT_FLOW (NUMODPAIR) =TRAFF I C 
END  IF 
END  DO 
ELSE 

WRITE (6, *) 'ORIGIN  IS  OUT  OF  BOUNDS,  MAX  ORIGIN=’,NN 
GO  TO  300 
END  IF 

STARTOD ( I ) = NUMODP AIR-NOD+ 1 
END  DO 
C 

C  OUTPUT  OF  CONNECTIVITY  DATA  FOR  DIRECT  INPUT  INTO  'MULTI FLO' 

C  COMMON  BLOCKS 

C 

WRITE (1,  *)NN 
DO  1=1, NN 

WRITE (1,*) FRSTOU (I) , LASTOU (I) 

END  DC 

WRITE  (1 ,  *)  NA 
DO  1=1, NA 

WRITE (1, «) STARTNODE (I) , ENDNODE ( I ) , BITRATE (I ) 

END  DO 


noon 


OUTPUT  OF  OD  TRAFFIC  DATA  FOR  DIRECT  INPUT  INTO  MULTIFLO 
COMMON  BLOCKS 

WRI TE ( 2 , * ) NUMCOMMOD 
DO  1=1, NUMCOMMOD 

WRITE (2, * ) ORGID (I) ,STARTOD(I) 

END  DO 

WRITE (2, * ) NUMODPAIR 
DO  1=1, NUMODPAIR 

WRITE (2  *) DEST (I) , INPUT_FLOW (I) 

END  DO 

STOP 

END 


APPENDIX  II:  MIJLTIFL01  Code 


The  only  differences  between  MULTIFLO  and  MIILTIFL01  are  in  the 
DRIVER  program  and  in  the  main  algorithm  subroutine  MIJLTIFLO.  These  two 
routines  called  DRIVER1  and  MULTIFL01 .  are  listed  below. 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C  DRIVER1 

c 

C  'DRIVER1  *  IS  A  SIMPLE  EXECUTIVE  TO  INVOKE  THE  ‘MULTIFLOl*  COMMODITY 

C  ROUTING  PROGRAM.  ' DRIVER1 *  INVOKES  SUBPROGRAM  'LOAD'  TO  READ 

C  DATA  INTO  'MULTIFLOl'  INPUT  COMMON  BLOCKS.  FILES  READ  BY 

C  'LOAD'  ARE  CREATED  BY  A  TERMINAL  SESSION  WITH  THE  USER  FOR 

C  NETWORK  DEFINITION  THROUGH  THE  USE  OF  PROGRAM  'SETUP'. 

C 

C  EXECUTION  STEPS  FOR  PROGRAM  'DRIVER1' 

C 

C  1)  ASSIGN  FORTRAN  UNIT  01  AS  CREATED  BY  PROGRAM  'LOAD' 

C  2)  ASSIGN  FORTRAN  UNIT  02  AS  CREATED  BY  PROGRAM  'LOAD' 

C  3)  ASSIGN  FORTRAN  UNIT  06  AS  A  DESIGNATED  OUTPUT  FILE 

C 

C  E .G. : 

C  $  ASSIGN  NETWORK . DAT  FOROOl 

C  $  ASSIGN  TRAFFIC.DAT  FOR002 

C  $  ASSIGN  OUTPUT . DAT  FOR006 

C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

PROGRAM  DRIVER1 
C 

C  LOAD  FORTRAN  UNIT  01  AND  FORTRAN  UNIT  02  FROM  DISK  AS  CREATED 

C  FROM  PROGRAM  'SETUP' 

C 

INCLUDE  'PARAM.DIM' 

INCLUDE  ' PATHS . BLK ' 

INCLUDE  'NETWRK.PRM' 

INCLUDE  ' CONVRG . PRM ' 

INTEGER  COMMOD ITY.ORIGIN.DE STOD , OD , PATH 
CALL  LOAD 
C 

C  EXECUTE  THE  'MULTIFLOl'  NETWORK  ALGORITHM.  'MULTIFLOl'  SCHEDULES 

C  ITS  OWN  OUTPUTS  TO  FORTRAN  UNIT  06  ON  EACH  ITERATION 

C 

C  INITIALIZE  THE  TIMER 

CALL  LIB$INIT_TIMER 
CALL  MULTIFLOl 

C  RECORD  THE  COMPUTATION  TIME 

CALL  LI B$SHOW_TIMER 
C 

C  PRINT  MAX  LINK  UTILIZATION  (RELEVANT  FOR  M/M/1  QUEUEING  DELAY 

C  OPTIMIZATION) 

C 

UMAX=0 .0 
DO  100  1=1. NA 

UMAX=MAX (UMAX, FA (I) /BITRATE (I) ) 

100  CONTINUE 

WRITE (6 , * ) 'MAXIMUM  LINK  UTILIZATION' 

WRITE (6,*) UMAX 

C  - 

C  PRINT  FINAL  PATH  FLOW  INFO 

C 

WRITE (6,*) 'ORIGIN  /  DESTINATION  /  PATH  #  /  PATH  FLOW’ 

DO  1000  COMMOD I TY= 1 , NUMCOMMOD 
ORIGIN=ORCID (COMMODITY) 

DO  500  OD=STARTOD (COMMODITY) , STARTOD (COMMODITY* 1) -1 


1 


500 

1000 


DESTOD-DEST (OD) 

PATH=OD 

DO  WHILE  (PATH.OT.O) 

WRITE (6, *)ORIGIN,DESTOD,PATH,FP (PATH) 
PATH=NEXTPATH (PATH) 

END  DO 
CONTINUE 
CONTINUE 
STOP 
END 


n  n  n 


3 


MULTICOMMODITY  FLOW  ALGORITHM  BASED  ON  A  PATH  FLOW  FORMULATION 
UPDATES  THE  PATH  FLOWS  OF  OD  PAIRS  ONE  AT  A  TIME  ACCORDING  TO 
AN  ITERATION  OF  THE  PROJECTION  TYPE. 

DEVELOPED  BY  DIMITRI  BERTSEKAS,  BOB  GENDRON,  AND  WEI  K  TSAI 
BASED  ON  THE  PAPERS: 


1) 


2) 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  MULTIFLOl 

c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 
c 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

SUBROUTINE  MULTIFLOl 

c 

IMPLICIT  NONE 


3) 


4) 


BERTSEKAS, D. P. ,  *'A  CLASS  OF  OPTIMAL  ROUTING  ALGORITHMS 
FOR  COMMUNICATION  NETWORKS".  PROC.  OF  5TH  I TERNAT I ONAL 
CONFERENCE  ON  COMPUTER  COMMUNICATION  (ICCC-80) , 

ATLANTA,  GA.,  OCT.  1980,  PP. 71-76. 

BERTSEKAS, D.P.  AND  GAFNI , E .M. ,  "PROJECTION  METHODS 
FOR  VARIATIONAL  INEQUALITIES  WITH  APPLICATION  TO 
THE  TRAFFIC  ASSIGNMENT  PROBLEM",  MATH.  PROGR.  STUDY. 17, 
D. C. SORENSEN  AND  J.-B.  WETS  (EDS),  NORTH-HOLLAND, 
AMSTERDAM, 1982,  PP.  139-159. 

BERTSEKAS, D. P. ,  "OPTIMAL  ROUTING  AND  FLOW  CONTROL 
METHODS  FOR  COMMUNICATION  NETWORKS",  IN  ANALYSIS  AND 
OPTIMIZATION  OF  SYSTEMS,  (PROC.  OF  5TH  INTERNATIONAL 
CONFERENCE  ON  ANALYSIS  AND  OPTIMIZATION,  VERSAILLES, 
FRANCE),  A.  BENSOUSSAN  AND  J.L.  LIONS  (EDS), 
SPRINGER-VERLAG,  BERLIN  &  NY, 1982,  PP.  615-643. 

BERTSEKAS, D.P.  AND  GAFNI ,  E.M.,  "PROJECTED  NEWTON 
METHODS  AND  OPTIMIZATION  OF  MULTICOMMODITY  FLOWS", 

IEEE  TRANSACTIONS  ON  AUTOMATIC  CONTROL,  DEC.  1983. 


*************** 


INCLUDE  COMMON  BLOCKS 


**************************** 


INCLUDE  'PARAM.DIM' 

INCLUDE  'NETWRK.PRM' 

INCLUDE  ' CONVRG . PRM ' 

INCLUDE  ' PATHS. BLK’ 

C 

C  NODE  ARRAYS  (LENGTH  NN)  : 

C  FRSTOU (NODE)  -  FIRST  ARC  OUT  OF  NODE 

C  LASTOU (NODE)  -  LAST  ARC  OUT  OF  NODE 

C  NOTE:  THE  ARC  LIST  MUST  BE  ORDERED  IN  SEQUENCE  SO 

C  THAT  ALL  ARCS  OUT  OF  ANY  NODE  ARE  GROUPED  TOGETHER 

C 

C  ARC  ARRAYS  (LENGTH  NA)  : 

C 

C  FA  (ARC)  -  THE  TOTAL  FLOW  OF  ARC 

C  STARTNODE (ARC)  -  THE  HEAD  NODE  OF  ARC 

C  ENDNODE (ARC)  -  THE  TAIL  NODE  OF  ARC 

C 


o  n 


COMMODITY  LENGTH  ARRAYS  (LENGTH  NUMCOMMOD) : 


C 

C 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


/ 


1 


c 

c 

c 


c 


c 

c 


c 

c 

c 


c 


c 

r 


ORGID (COMMODITY)  -  THE  NODE  ID  OF  THE  ORIGIN  OF  COMMODITY 
STARTOD (COMMODITY)  -  THE  STARTING  OD  PAIR  IN  THE  ODPAIR  LIST 
CORRESPONDING  TO  THE  ORIGIN  IN  POSITION  RANK 
NOTE:  THIS  SCHEME  ASSUMES  THAT  OD  PAIRS  ARE  LISTED  IN  SEQUENCE 
I.E.  THE  OD  PAIRS  CORRESPONDING  TO  THE  COMMODITY  ONE 
ARE  LISTED  FIRST.  THEY  ARE 

FOLLOWED  BY  THE  OD  PAIRS  OF  THE  COMMODITY  TWO 
AND  SO  ON. 

ODPAIR  ARRAYS  (LENGTH  NUMOD) : 

DEST (OD)  -  GIVES  THE  DESTINATION  OF  ODPAIR  OD 
INPUT_FLOW (OD)  -  GIVES  THE  INPUT  TRAFFIC  OF  ODPAIR  OD 

PATH  ARRAYS  (LENGTH  DYNAMICALLY  UPDATED) : 

PATHID (PATH)  -  POINTER  TO  THE  BLOCK  DESCRIBING  PATH 
IN  THE  PATH  DESCRIPTION  ARRAY 

NEXTPATH (PATH)  -  THE  NEXT  PATH  FOR  THE  SAME  OD  PAIR  FOLLOWING 
PATH.  IT  EQUALS  0  IF  PATH  IS  THE  LAST  FOR  THAT  OD  PAIR 
FP  (PATH)  -  THE  FLOW  CARRIED  BY  PATH 

PATH  DESCRIPTION  LIST  ARRAY  (LENGTH  DYNAMICALLY  UPDATED) 

PDESCR (INDEX)  -  THIS  LONG  ARRAY  EXPLICITLY  DESCRIBES  ALL 
ACTIVE  PATHS.  FOR  ANY  PATH,  PATHID (PATH)  IS  A  POINTER 
TO  PDESCR.  IT  GIVES  THE  ELEMENT 

OF  THE  PDESCR  ARRAY  CONTAINING  THE  #  OF  ARCS  IN  THE  PATH 
(CALL  IT  NUMARC) .  THE  ELEMENTS  PATHID (PATH) -NUMARC  TO 
PATHID  (PATH) -1  OF  THE  ARRAY  PDESCR  CONTAIN  THE  ARCS  THAT 
MAKE  UP  PATH  STARTING  FROM  THE  DESTINATION  AND  GOING  TOWARDS 
THE  ORIGIN  OF  PATH. 

***************  LOCAL  VARIABLE  DEFINITIONS  ************************ 

INTEGER *2  PDESCR (NINDEX) 

PATH  DESCRIPTION  ARRAY  -  CONTAINS  EXPLICIT 
DESCRIPTION  OF  ALL  ACTIVE  PATHS. 

LOGICAL  SPNEW 

LOGICAL  INDICATING  A  NEW  PATH  FOUND 
LOGICAL  SAME 

LOGICAL  INDICATING  A  NEW  SHORTEST  PATH  ALREADY  EXISTING 
INTEGER  NODE 

NODE  IDENTIFIER 
INTEGER  DESTOD 

THE  DESTINATION  NODE  OF  AN  OD  PAIR 
INTEGER  ARC 

DO  LOOP  INDEX  FOR  ARCS  -- 

INTEGER  PATH 

A  PATH  INDEX 
INTEGER  NUMLIST 

TOTAL  NUMBER  OF  ACTIVE  PATHS  FOR  OD  PAIR  UNDER  CONSIDERATION 
INTEGER  ITER 

SPECIFIC  ITERATION 
INTEGER  N1,N2 

TEMPORARY  VARIABLES 
REAL  MINFDER 

THE  LENGTH  FOR  A  SHORTEST  PATH 
REAL  MINSDER 

THE  SECOND  DERIVATIVE  LENGTH  FOR  THE  SHORTEST  PATH 
REAL  TMINSDER 

TEMPORARY  VALUE  FOR  SECOND  DERIVATIVE  LENGTH  OF  SHORTEST  PATH 


o  u>  nooonnno 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


REAL  I NCR 

TOTAL,  SHIFT  OF  FLOW  TO  THE  MINIMUM  FIRST  DERIVATIVE  LENGTH  PATH 
REAL  PATH I NCR 

SHIFT  OF  FLOW  FOR  A  GIVEN  PATH 
REAL  FLOW 

FLOW  FOR  A  PATH 
REAL  FDER 

THE  ACCRUED  LENGTH  ALONG  A  PATH 
REAL  SDER 

THE  ACCRUED  SECOND  DERIVATIVE  LENGTH  ALONG  A  PATH 
REAL  TEMP  ERROR 

TEMPORARY  STORAGE  FOR  CONVERGENCE  ERROR 
REAL  FDLENGTH (NMAX ITER) 

ARRAY  OF  LENGTHS  OF  PATHS  FOR  AN  OD  PAIR 
REAL  SDLENGTH (NMAX I TER) 

ARRAY  OF  SECOND  DERIVATIVE  LENGTHS  OF  PATHS  FOR  AN  OD  PAIR 
INTEGER  PATHLIST (NMAXITER) 

ARRAY  OF  ACTIVE  PATHS  FOR  AN  OD  PAIR 
INTEGER  COMMODITY 

DO  LOOP  INDEX  FOR  THE  OD  PAIR  ORIGINS 
INTEGER  ORIGIN 

SPECIFIC  ORIGIN 
INTEGER  I 

DO  LOOP  INDEX 
INTEGER  OD 

OD  DO  LOOP  INDEX 
INTEGER  K 

DO  LOOP  INDEX 
INTEGER  SHORTEST 

THE  SHORTEST  PATH 
INTEGER  INDEX 

THE  CURRENT  LAST  ELEMENT  OF  THE  ARRAY  PDESCR 
INTEGER  POINT 

POINTER  TO  PDESCR 
INTEGER  NUMARC 

#  OF  ARCS  IN  A  PATH 
LOGICAL  MEMBER (NNA) 

LOGICAL  FOR  AN  ARC  INCLUDED  IN  THE  SHORTEST  PATH 
REAL  DLENGTH 

DIFFERENCE  IN  PATH  LENGTHS  FOR  THE  TRAFFIC 
REAL  D1CAL 

ARC  LENGTH 
REAL  D2CAL 

DERIVATIVE  OF  ARC  LENGTH 


**********************  EXECUTABLE  CODE 


***************************** 


*  INITIALIZATION 

***************************************** 


DO  5  ARC=1 ,NA 
FA  (ARC)  =0.0 
CONTINUE 

DO  1=1 , NUMODPAIR 

FP (I) =INPUT_FLOW (I) 

ENDDO 

STARTOD (NUMCOMMOD+ 1 ) =NUMODPAIR+ 1 
NUMPATH=0 


INDEX=0 

NUMITER=1 

DO  100  COMMODI TY=1 , NUMCOMMOD 
ORIGIN=ORGID (COMMODITY) 

CALL  SP (ORIGIN, COMMODITY) 

C 

C  LOOP  OVER  OD  PAIRS  OF  COMMODITY 

C 

Nl=STARTOD (COMMODITY) 

N2=STARTOD (COMMODITY+1)  -I 
DO  50  OD-Nl ,N2 

NUMPATH=NUMP ATH+ 1 
NEXTPATH (NUMPATH) =0 
FLOW=FP (NUMPATH) 

INDEX=INDEX+1 
NUMARC=0 
NODE=DEbT (OD) 

DO  WHILE  (NODE. NE . ORIGIN) 

ARC=PA(NODE) 

FA (ARC) =FA (ARC) +FLOW 
PDESCR (INDEX) =ARC 
NUMARC=NUMARC+ 1 
INDEX=INDEX+1 
NODE=STARTNODE (ARC) 

END  DO 

PATHID  (NUMPATH)  =  INDEX 
PDESCR  (INDEX)  =NUMARC 
50  CONTINUE 

100  CONTINUE 

C 

C  INITIALIZE  MEMBER  ARRAY 

C 

DO  70  ARC=1,NA 

MEMBER (ARC) = . FALSE . 

70  CONTINUE 

C 

C  INITIALIZE  THE  TOTAL  DELAY 

C 

CALL  DELAY (DTOT (NUM1TER) ) 

C 

C  OUTPUT  THE  CURRENT  INFORMATION  TO  DISK 

C 

CALL  PRFLOW 
C 

C  ********************************************** 

C  *  END  OF  INITIALIZATION 

C  ********************************************** 

C 

C  *****  START  NEW  ITERATION  ***** 

C 

110  NUMITER=NUMITER+1 

CURERROR=0 
C 

C  ****  LOOP  OVER  ALL  COMMODITIES  **** 

C 

DO  1000  COMMODI TY=1 .NUMCOMMOD 
OR I G I N=ORGI D (COMMODITY) 

CALL  SP  (OR  I  GIN,  COMMODITY) 

C 

r  LOOP  OVER  OD  PAIRS  OF  COMMODITY 


c 


N 1 = STARTOD ( COMMOD I TY ) 

N  2=STARTOD  ( COMMOD I  TY  + 1 )  - 1 
DO  500  OD=Nl , N2 
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C 

C  CHECK  IF  THERE  IS  ONLY  ONE  ACTIVE  PATH  AND  IF  SO  SKIP 

C  THE  ITERATION 

C 

IF  (NEXTPATH(OD) .EQ.O)  THEN 
NODE=DEST (OD) 

POINT=PATHID (OD) 

NUMARC=PDESCR (POINT) 

DO  150  I =POINT-NUMARC , POINT- 1 
ARC=PDESCR ( I ) 

IF  (ARC. NE .PA (NODE) )  GO  TO  180 
NODE=STARTNODE (ARC) 

150  CONTINUE 

GO  TO  500 
END  IF 


C 

180  CONTINUE 

C 

C  MARK  THE  ARCS  OF  THE  SHORTEST  PATH 

C 

DESTOD=DEST (OD) 

NODE=DESTOD 

DO  WHILE  (NODE. NE. ORIGIN) 

ARC=PA (NODE) 

MEMBER (ARC) = . TRUE . 

NODE=STARTNODE (ARC) 

END  DO 
C 
C 

C  GENERATE  LIST  OF  ACTIVE  PATHS  FOR  OD  PAIR 

C 

NUMLI ST=1 
PATHLIST(l)=OD 
PATH=NEXTPATH (OD) 

DO  WHILE  (PATH.GT.O) 

NUMLI ST=NUMLIST+1 
PATHLI ST (NUMLI ST) =PATH 
PATH=NEXTPATH (PATH) 

END  DO 
C 

C  DETERMINE  1ST  &  2ND  DERIVATIVE  LENGTH  OF  ACTIVE  PATHS 

C  ALSO  DETERMINE  WHETHER  THE  CALCULATED  SHORTEST  PATH 

C  IS  ALREADY  IN  THE  LIST 

C 

SPNEW= . TRUE . 

DO  200  K=l. NUMLI ST 
SAME= .TRUE . 

FDER=0 

SDER=0 

TMINSDER=C 


PATH=PATHLIST (K) 
POINT=PATHID (PATH) 
NUMARC=PDESCR (POINT) 


DO  210  I=POINT-NUMARC, POINT-1 
ARC=PDESCR ( I ) 

CALL  DERI VS (COMMODITY, FA (ARC) , ARC , D1CAL , D2CAL) 


FDER=FDER+D1CAL 
IF  (. NOT. MEMBER (ARC) )  THEN 
SDER=SDER+D2CAL 
SAME=. FALSE. 

ELSE 

SDER=SDER-D2CAL 
TMINSDER=TMINSDER+D2CAL 
END  IF 
CONTINUE 
IF  (SAME)  THEN 

SPNEW=. FALSE. 

SHORTEST=PATH 
FDLENGTH (K) =FDER 
MINFDER=FDER 
MINSDER=TMINSDER 

ELSE 

FDLENGTH (K) =FDER 
SDLENGTH (K) =SDER 
END  IF 
CONTINUE 

***  INSERT  SHORTEST  PATH  IN  PATH  LIST  IF  IT  IS  NEW  *** 

IF  (SPNEW)  THEN 

NUMPATH^NUMPAIH  + 1 
SHORTEST=NUMPATH 

NEXTPATH (PATHLI ST (NUMLI ST) ) =NUMPATH 

NEXTPATH (NUMPATH) =0 

MINFDER=0 

MINSDER=0 

INDEX=INDEX+1 

NUMARC=0 

NODE=DESTOD 

DO  WHILE  (NODE. NE .ORIGIN) 

ARC=PA(NODE) 

PDE SCR  (INDEX)  =ARC 
NUMARC=NUMARC+ 1 
INDEX=INDEX+1 

CALL  DERI VS (COMMODITY, FA (ARC) , ARC,D1CAL,D2CAL) 

MINFDER=MINFDER+D1CAL 

MINSDER=MINSDER+D2CAL 

NODE=STARTNODE (ARC) 

END  DO 

PATHID (NUMPATH) =INDEX 
PDESCR (INDEX) =NUMARC 
END  IF 

****  UPDATE  PATH  &  LINK  FLOWS  **** 


INCR=0 

TEMP ERROR =0 

DO  250  K=1 .NUMLIST 

DLENGTH=FDLENGTH  (K)  -MINFDER 
IF  (DLENGTH.GT.O)  THEN 
PATH=PATHLIST(R) 

FLOW=FP  (PATH) 

IF  ( (FLOW. EQ. 0.0) .AND. (K.GT.l))  THEN 
NEXTPATH  (PATHLI  ST  (K- 1 )  )  =  NEXTPATH  (PATH) 
GO  TO  250 
END  IF 


PATHINCR=DLENGTH/ (SDLENGTH (K) +MINSDER) 

IF  (FLOW . LE . PATHINCR)  THEN 
FP (PATH) =0.0 
PATHI NCR=FLOW 
ELSE 

FP (PATH) =FLOW-PATHINCR 
END  IF 

I NCR= I NCR +P ATHI NCR 

TEMP  ERROR=TEMP ERROR + FLOW  *  DLENGTH/FDLENGTH (K) 
POINT=PATHID (PATH) 

NUMARC=PDESCR (POINT) 

DO  220  I =POI NT-NUMARC , POINT- 1 
ARC=PDESCR (I) 

FA (ARC) =FA (ARC) -PATHINCR 
220  CONTINUE 

END  IF 

250  CONTINUE 

C 

C 

C  ***  UPDATE  THE  ERROR  CRITERION  *** 

C 

CUR ERROR =AMAX1 (CURERROR , TEMPERROR/INPUT_FLOW (OD) ) 
C 

c  ****  UPDATE  FLOWS  FOR  SHORTEST  PATH  **** 

C 

FP  (SHORTEST) =FP (SHORTEST) +INCR 
POINT=PATHID (SHORTEST) 

NUMARC=PDESCR (POINT) 

DO  300  I =POI NT-NUMARC , POINT- 1 
ARC=PDESCR (I) 

FA (ARC) =FA (ARC) + I NCR 
MEMBER (ARC) = . FALSE . 

300  CONTINUE 

C 

500  CONTINUE 

C 

C  *****  end  OF  LOOP  FOR  OD  PAIRS  CORRESPONDING  TO  COMMODITY 

c  *****  UPDATE  TOTAL  DELAY 

C 

CALL  DELAY (DTOT (NUMITER) ) 

C 

1000  CONTINUE 
C 

C  CHECK  IF  THE  #  OF  ACTIVE  PATHS  EXCEED  THE  ALLOCATED  NUMBER 

C 

IF  (NUMPATH . GT . NNUMPATH)  THEN 

WRITE (6, *) ’MAX  #  OF  ALLOCATED  PATHS  EXCEEDED' 

STOP 
END  IF 

IF  (INDEX. GT.NINDEX)  THEN 

WRITE (6,*) 'DIMENSION  OF  PDESCR  ARRAY  EXCEEDED' 

STOP 
END  IF 
C 

C  OUTPUT  THE  CURRENT  SOLUTION  TO  DISK 

C 

CALL  PRFLOW 
C 

C  *****  end  OF  ITERATION  ***** 


nnnn 


***  IF  THE  ERROR  IS  SMALLER  THAN  TOL,  OR  THE  LIMIT  ON 
THE  NUMBER  OF  ITERATIONS  IS  REACHED  RETURN 
ELSE  GO  FOR  ANOTHER  ITERATION 

IF  ( (CUR ERROR .LT. TOL) .OR. (NUMITER.EQ.MAXITER) )  THEN 

WRITE (6, *) ’FINAL  STORAGE  OF  PATH  DESCRIPTION  LIST' 
WRITE  (6,  *)  INDEX 
RETURN 

ELSE 

GO  TO  110 
END  IF 
C 

END 

C  **************  END  OF  MULTIFLOl  **************** 


